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ABSTRACT 

The dynamics of magnetic fields in closed regions of solar and stellar coronae are investigated with a 
reduced magnetohydrodynamic (MHD) model in the framework of Parker scenario for coronal heating. 
A novel analysis of reduced MHD equilibria shows that their magnetic fields have an asymmetric 
structure in the axial direction with variation length-scale zi ^ iBo/h^ where Bq is the intensity 
of the strong axial guide field, b that of the orthogonal magnetic field component, and £ the scale 
of b. Equilibria are then quasi-invariant along the axial direction for variation scales larger than 
approximatively the loop length Z£ ^ Lz^ and increasingly more asymmetric for smaller variation scales 
^ Bz. The critical length zg ^ Lz corresponds to the magnetic field intensity threshold b ^ £Bq/Lz. 
Magnetic fields stressed by photospheric motions cannot develop strong axial asymmetries. Therefore 
fields with intensities below such threshold evolve quasi-statically, readjusting to a nearby equilibrium, 
without developing nonlinear dynamics nor dissipating energy. But stronger fields cannot access 
their corresponding asymmetric equilibria, hence they are out-of-equilibrium and develop nonlinear 
dynamics. The subsequent formation of current sheets and energy dissipation is necessary for the 
magnetic field to relax to equilibrium, since dynamically accessible equilibria have variation scales 
larger than the loop length zi > Lz^ with intensities smaller than the threshold b < £Bq/Lz. The 
dynamical implications for magnetic fields of interest to solar and stellar coronae are investigated 
numerically and the impact on coronal physics discussed. 

Keywords: magnetohydrodynamics (MHD) — stars: activity — stars: solar-type — Sun: corona — 
Sun: magnetic topology — turbulence 


1. INTRODUCTION 

Solar observations show a close association between 
magnetic field strength and coronal activity. In combi¬ 
nation with the ability of photospheric motions to stress 
the field, these are the two key elements to understand 
the observed coronal X-ray activity of the Sun, of all 
late-type main sequence stars, and more in general of 
stars wit h a magnetized cor ona and an outer convective 
envelope (jGudelll2QQ4 l2QQ^ . 

It has long been proposed that the work done by 
photospheric motions on magnetic field line footpoints 
can transform mechanical energy into magnetic energy 
and transfer it in the upper corona. The photospheric 
(horizontal) velocity can be split into irrotational and 
solenoidal components. Only its solenoidal {incompress¬ 
ible) component has a non-vanishing vorticity and can 
then twist the magnetic field lines, injecting magnetic 
e nergy into the corona . 

I Gold fc Hovld (|196Ql ) conjectured that the magnetic 
field would proceed through a sequence of force-free equi¬ 
libria while photospheric vortices twist the field lines, 
and the stored energy could subsequently be released 
when two flux tubes with similarly twisted field lines 
come into contact with each other, or when the configu¬ 
ration would become so mehow uns t able through an un¬ 
determ ined mechanism (|Goldlll9^ . I St nr rock fc Uchidl] 
compute the energy flux into the corona due to 
the work done by random photospheric vortical motions 
on the magnetic field. They find that the correlation 
time of photospheric motions must be of the order of the 
observed photospheric timescales (5-8 minutes) or longer 


to obtain an energy flux large enough to sustain an ac¬ 
tive corona, otherwise for shorter correlation times the 
resulting twisted field is too small. But the magnetic en¬ 
ergy is still supposed to be stored in a force-free field in 
equilibrium, and no physical mechanism able to dissipate 
this energy and heat the corona is envisioned. 

Energy stored in a magnetic field in equilibrium, that 
subsequently becomes unstable and r eleases its energy, is 
the common thread of fl are models (jShibata fc Magarl] 
1201 ll: iMartin et al.|[2QT^ . with the processes leading to 
the pre-flare magnetic energy storage and its subsequent 
fast release strongly debated. On the other hand this 
picture does not appear apt to describe the dynamics of 
the long-lived slender X-ray bright loops, that in com¬ 
parison to a flare show little dynamics from their large- 
scales down to the smallest resolved scale (^ 150 km) 
of current state-of-the-art X-r ay and EUV ima gers on 
board Hinode, SDO and Hi-C (|Peter et al.ll2013[ ). While 
the pre-flare magnetic structure is destroyed during the 
flare, the large-scale magnetic topology of the loops where 
the basic heating occurs, and that make the corona shine 
steadily in X-ray, is not strongly modified on compara¬ 
ble timescales. This suggests that the energy deposition 
must occur at very small scales yet observationally un¬ 
resolved. Eurthermore the energy reservoir that supplies 
dissipation should consist of magnetic field fluctuations 
(with vanishing time-average, but non vanishing time- 
averaged rms) that adds up to the strong axial magnetic 
fi eld that defines the l oop. 

iParkerl (|1972l . 1 19881 . 11994 l2012f ) was the first to sug¬ 
gest that, in contrast to previous quasi-static models, 
the magnetic field brought about by photospheric vorti- 
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cal motions would be in dynamical non-equilibrium in 
the case of interest to coronal heating. Furthermore 
the relaxation of this interlaced fields toward equilibrium 
would necessarily involve the formation of current sheets. 
The energy dissipation would then occur at small scales 
in the fashion of small impulsive heating events, so-called 
nanoflares^ a picture broadly used for the thermodynam¬ 
ical m odeling of the closed corona (jKlimchukll2QQ6l: iReald 

IMl . 


Using a simplified Cartesia n model with a stron g guide 
field threading a coronal loop. lParkerl(jl972L 1 19791) argued 
at first that a magnetic field could be in equilibrium only 
if it were invariant along 2 ; (the axial direction). Due 
to the complex and disordered nature of photospheric 
motions the induced interlaced magnetic field would not 
be invariant, and therefore not in equilibrium. Next, 
counterexamples of magnetostatic e quilibria that are not 
invariant along ^ were p rovided (|Rosner fc Knoblochl 
119821: lBogovavlenskiill2QQQl), and analytica l investigations 
(Ivan Ballegooiienl 119851: lAntioch^ 119871: iCowlev et al.l 
II 997 D argued that smooth photospheric motions cannot 
lead to the formation of current sheets, whereas only a 
discontinuous velocity field can form discontinuities in 
the coronal ma gnetic field. _ 

In particular Ivan Ballegooiienl (|1985D showed that the 
equilibria are the solutions of the two-dimensional (2D) 
Euler equation that in general are not ^-invariant, thus 
inferring that the field would evolve quasi-statically, con¬ 
tinuously readjusting to a nearby force-free equilibrium 
without developing nonlinear dynamics no r forming cur¬ 
rent sheets. Reac hing opposite conclusions iParkerl ()l988L 
I 1994 I2QQQII 2 QI 2 D pointed out that almost all field line 
topologies relevant to the solar corona have a different 
structure from the solutions of the Euler equation, so 
that the magnetic field would be still in dynamical non- 
equilibrium. 

Reduced magnetohydrodynamics (MHD) numerical 
simulations with a continuous smooth velocity forcing at 
the boundaries show that the dynamics can be seen as a 
particula r instance of magnetically dominated MHD tur¬ 
bulence (iDmitruk fc Gon^ 119991: iDmitruk et al.l I2QQ3I: 


iRappazzo et al.l l2QQ7l 


models ( Einaudi et al 


2QQ8D as proposed by ear l ier 20 
119961: IDmitruk fc G6me3 I1997D . 


suggesting that in the forced case the magnetic field 
is in dynamical non-equilibrium rather than close to a 
quasi-static evolution. Similar dynamics are also dis¬ 


played b y boundary driven simulation s in the cold plasma 
regime (jHendrix fc van Hovenl 1 19961) and in the fully 
compressible MHD c ase (jGalsgaard fc Nordlundl 119961: 
iDahlburg et al.l I2Q12I ) . Eurthermore IRappazzo fc Vellil 
dmH) have shown that while velocity fluctuations are 
much smaller than magnetic fluctuations, spectral en¬ 
ergy fluxes toward smaller scales are akin to those of 
a standard cascade with magnetic and kinetic ener¬ 
gies in equipartition, except for kinetic energy fluxes 
that are negligible. This implies that at scales smaller 
than those directly shuffled by photospheric motions, the 
small velocity field is created and shaped by the unbal¬ 
anced Lorentz force of the out-of-equilibrium magnetic 
fields that in turn creates small scales in the magnetic 
field by distorting magnetic i slands and push i ng fiel d 
lines together. Add itionally iGeorgoulis et al.l (jl998D : 
IDmitruk et al ] (IT^ have established a link between 
boundary driven simulations and observed statistics of 


coronal activity. Indeed the bursts in dissipation dis¬ 
played by the system, that correspond to the formation 
and dissipation of current sheets, follow a power law be¬ 
havior in total energy, peak dissipation and duration with 
indexes not far from those determined observationally in 
X-rays. 

Recently IWilmot-Smith et al.l (|2QQ9[ ) have shown that 
the relaxation of a slightly braided magnetic field (“pig¬ 
tail” braid) appears to evolve quasi-statically, with no 
formation of current sheets, toward an equilibrium where 
only large-scale current layers of thickness much larger 
than the resolution scale are observed. This result would 
seem in contrast with Parker’s hypothesis, the results of 
the forced numerical simulations discussed in the previ¬ 
ous paragraph, and the recent results supporting the de- 
velopm ent of finite time singularities in the cold plasma 
regime (ILowll 2 ni.llM^ . 

To get fu rther insight on the dynam ics of coronal mag¬ 
netic fields, IRappazzo fc Parl^ (|2Q13f ) analyzed reduced 
MHD numerical simulations of the relaxation of initial 
magnetic fields invariant along z and with different av¬ 
erage twists. They identified a critical intensity threshold 
for the magnetic field. This is explained heuristically as 
due to a balance between different field line tension forces 
for weak fields, while such a balance cannot be attained 
by stronger fields. The non-equilibrium of stronger fields 
stems from this force unbalance, and drives the relax¬ 
ation forming current sheets and dissipating energy. On 
the contrary weaker fields show little dynamics with no 
energy dissipation, confirming that they are essentially in 
equilibri um. Such threshold can expl ain qualitatively the 
result of IWilmot-Smith et al.l (I2QQ9D . although a quan¬ 
titative comparison cannot be made because the inte¬ 
grated equations (magneto-frictional relaxation vs. re¬ 
duced MHD) and initial topologies differ. 

The magnetic i ntensi ty threshold found by 
IRappazzo fc Parl^ (|2Q13D implies that a critical 
twist exists above which dynamics develop, and below 
which t he sy stem remains very close to equilibrium. 
IParkerl (| 19881) had conjectured that a critical twist is 
necessary to explain t he observationally inferre d energy 
flux in active regions (|Withbroe fc No\^ll977l) . In fact 
the energy flux injected in the corona by photospheric 
motions is t he a verage Poynting flux {Sz) = Bq (uph • b) 
(see Section na Equation [7]) that depends not only on 
the photospheric velocity Uph and the axial guide field 
Ho 5 but also on the dynamic magnetic field component 
b, with stronger intensities corresponding to higher 
average twists. Thus if nonlinear dynamics were to 
develop for weak field intensities, energy dissipation 
would keep too low the value of b, and consequently 
the flux {Sz)‘ This argument is further reinforced 
by the fact that current sheets thickness decreases at 
least exponentially in time when nonlinear dy namics 
develo p, reaching t he Sweet-Parker thickness (|Sweetl 
[T^ IParkerl on id eal timescales (about on e 

Alfven crossing time IRappazzo fc Parl^ I2Q13D . 
that current sheets are unstable to tearing modes with 
“ideal” growth rates (i.e., of order I/t a) already at 
thicknesses larger than Sweet-Parker (iPucci fc Vellil 
|2Q14 iTenerani et al.ll2QT5l: iLandi et al.l 120151) . and that 
magnetic reconnection rates can b e very fast in plasmas 


with high Reynolds numbers ( 

Lazarian & Vishniac 

119991: iLoureiro et al.ll2007l: iLapenta 

I 2 OO 8 I: ILoureiro et al. 
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2009 : lUzdenskv et al.l 120101: iHuang fc Bhattacharied 
201C: iBeresnvakl I2013D and in th e collisionless regime 
(jShay et al.l 119991: iBirn et ^l2001f ). 

This paper is devoted to a more detailed discus- 
sio n and analysis of the nume rical simulations described 
by iRappazzo fc Parl^ (j2013[) , of additional simulations 
that extend our previous work to initial conditions non¬ 
invariant along z, and to a novel analysis of the struc¬ 
ture of the reduced MHD equilibria, with the goal to 
shed light on the topics outlined in this introduction and 
advance our understanding of coronal magnetic field dy¬ 
namics, their relationship to dynamic non-equilibrium, 
MHD turbulence, quasi-static evolution, current sheets 
formation and activity of solar and stellar coronae. 

The loop model along with initial and boundary condi¬ 
tions for the simulations are introduced in Section O The 
structure of the equilibria is analyzed in Section [H and 
the results of the numerical simulations are described in 
Section IH Finally results and conclusions are summa¬ 
rized in Section [5l together with a discussion of their 
impact on coronal physics. 

2. PHYSICAL MODEL 

A closed regi on of the solar corona is modeled, as in 
previous work (jRappazzo et al.ll2QQ7[ ). with a simplified 
cartesian geometry, uniform density p and a strong and 
homogeneous axial magnetie field Bq = Bq thread¬ 
ing the system. Plasmas in such configurations are 
well suited to be studied in the reduced MHD regime 
(|Zank fc Matthaeuslll992[ ). Introducing the velocity and 
magnetic field potentials (p and for which u = V(/:?x e^, 
b = X e^, vorticity uj = and the current den¬ 

sity j — — the nondimen sional reduced MHD equa¬ 
tions (jKadomtsev fc Pogutselll97^ lStranssl[l976[ ) are: 

dtip = [<p, Ip] + Bo + rjj^pip, ( 1 ) 

= [j, Ip] - [uj, <p] + Bo dzj + (2) 

The Poisson bracket of functions g and h is defined 
as [g,h] = dxgdyh - dygd^h (e.g., [j,^)] = b ■ Vj), 
and Laplacian operators have only orthogonal compo¬ 
nents. To render the equations nondimensional the 
magnetic field s are first expressed as Alfven velocities 
(6 ^ b/^/A'Kp)^ and then all velocities are normalized 
with 14* = 1 kms“^, a typical value for photospheric mo¬ 
tions. The domain spans 0 < x, ^ and ^ < z < Lz^ 

with = 1 and Lz = 10. Magnetic field lines are 
line-tied to a motionless photosphere at the top and bot¬ 
tom plates {z = 0 and 10), where a vanishing velocity 
u = 0 is in place. In the perpendicular (x-y) direc¬ 
tions a pseudo-spectral scheme with periodic boundary 
conditions and isotropic trun cation de-aliasing is used 
(2/3-rule. iCanuto et al.ll2006f ). while along 2 : a second- 
order finite difference scheme is implemented. The CFL 
(Courant-Friedrichs-Levy) condition is satisfied through 
an adaptive time-step. For a mor e detailed descrip t ion of 
the m odel and numerical code see IRappazzo et ahl (j2007L 
[200^1 . 

Dis sipative simulations use hyper-diffusion (|BiskamDl 
l2003f ). that effectively limits diffusion to the small scales, 

with n = 4 and i^n = hn = /Bn^ with 

corresponding to the Reynolds number for n = 1 (see 
IRappazzo et al.ir2QQ8f ). 


2.1. Initial and boundary eonditions 

Simulations are started at time t = 0 with a vanishing 
velocity u = 0 everywhere, and a uniform and homoge¬ 
neous guide field Bq. The orthogonal field b consists of 
large-seale Fourier modes, set expanding the magnetic 
potential in the following way: 


p,o=hoY.{2SJ 


1 CXrsm sin (k7-5772 ' X T 


kr 




(3) 


2ti 


ij ijm 

2tt 


with \^rsm = + s By) ^ -—mez, 

L±_ Lz 


and kr 


+ s^, 

L^ 


where the coefficients o^rsm and ^rsm are two indepen¬ 
dent sets of random numbers uniformly distributed be¬ 
tween 0 and 1. The orthogonal wave-numbers (r, 5 ) are 
always in the range 3 < (r^ -1- < 4, while the 

parallel amplitudes S'm (with = 1) are set to 

distribute the energy in different ways in the axial direc¬ 
tion. Given the orthogonality of the base used in Eq. m 
the normalization factors guarantee that for any choice 
of the amplitudes the rms of the magnetic field is set to 
^ — {^x = ^0, while for total magnetic energy 

Bm — ^ 0^/2 he., S’m is the fraction of magnetic 

energy in the parallel mode m. Two-dimensional (2D) 
configurations invariants along ^ are obtained consider¬ 
ing the single mode m = 0 with (^0 = 1. 


2.2. Energeties 

From equations with n = 1 and considering 

the kinetic and magnetic Reynolds numbers equal, the 
following energy equation can be obtained: 

I i(u2 + b2) = -V-(S + F)-l(y + ^2)^ (4) 

where S = Bx(uxB)is the Poynting vector, and 
F = (p -h u^/2)u — (cju + jh) X Bz/R is an orthogonal 
transport-related flux. Integrating equation (|1]) over the 
whole box the energy (E) equation is 



i.e., as expected, the global energy balance depends on 
the competition between the energy flowing into the com¬ 
putational box from the photospheric boundaries S and 
the ohmic and viscous dissipation. Because in the x- 
y planes periodic boundary conditions are implemented, 
and = 0, the only relevant component of the flux 
vectors is that of the Poynting vector along the axial di¬ 
rection Sz that, as B = Hoe;^ -h b, is given by 

5, = S ■ e, = -Bo (u ■ b). ( 6 ) 

Indicating the photospheric velocity fields at the top and 
bottom boundaries z = 0 and L with and for the 
integrated energy flux (i.e., the power) S we obtain 

S = Bq j da (u^ • b) — Ho J da (u^ • b) . (7) 

z=L z=0 


The injected energy power is proportional to Bq and de¬ 
pends on the dot product of the photospheric velocities 
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and the magnetic field b at the boundaries. But 
while and Bq have fixed values (in our simplified 
model), the magnetic field b is determined by the linear 
or nonlinear dynamics developing in the computational 
box. 

Since the magnetic field component b can often be 
considered quasi-invariant along ^ (as described in the 
following sections), as a shorthand we will indicate the 
difference between the boundary velocities with Wph = 
— u^, so that for quasi-invariant fields the Poynting 
flux can be approximated as {Sz) ^ SjP r\j BoUphb. 


3. EQUILIBRIA AND THEIR DYNAMIC 
ACCESSIBILITY: ANALYSIS 


As discussed in the introduction, the properties of the 
equilibri a of this system are pivotal to understand its dy¬ 
nami ses (|Parkerlll972L 1 1 988111994 Ivan Ballegooiien|[l985L 
1 1 9861) . therefore their structure is analyzed here in detail. 
It is shown that, depending on the ratio bo/Bo of the rms 
of the orthogonal component to the guide magnetic field 
intensity, the equilibria ean be approximately invariant 
along z or strongly asymmetrie. As shown in the follow¬ 
ing this can explain why fields with a twist below a eriti- 
eal value do not form strong current sheets, while they do 
at higher twists as conjectured bv iParken (|1988f ). Nev¬ 
ertheless unlike commonly thought, such equilibria are 
generally not linearly unstable for most conditions rele¬ 
vant to coronal loops, since they arise from a balance of 
forces in an asymmetrie and irregular topology. 

Neglecting velocity and diffusion terms, equilibria of 
Eqs. ([I])-© are given by B • Vj = 0. Since the total 
magnetic field B is given by B = BqSz + b(x, z) with 
h Bz = 0, the equilibrium condition can be written as: 


dz 



( 8 ) 


where the right-hand side term corresponds to the “2D 
perpendicular” Lorentz force component b • Vb, and the 
left hand side to the “parallel” B^dzh field line tension 
(the labels refer to their derivative, but both components 
are ort hogo nal to Rq, a more detailed discussion is in 
Section I01 prior to Equation [23]). 

Assigned b in an x-y plane, e.g., at the boundary z=0 
h{x^y^z = 0) = hi)(i{x^y)^ the integration of this equa¬ 
tion for 2 : > 0 allows to compute the corresponding equi¬ 
librium in the whole computational box 0 < ^ < 

Now consider the 2D Euler equation (jEulerlllTGlI ) 

| = -u.V„, (9) 

with V-u = 0. Introducing the velocity potential 0, then 
u(x,^,t) = V0(x,^,t) X ez^ and vorticity uj = —V‘^(j). 
The two equations © and © are identieal under the 
mapping 


t —}■ 2 :, 

< b (10) 

u , 


and consequently uj j jB^. 

The related 2D Navier-Stokes equation is obtained by 
adding to the right hand side of Equation © the dissi¬ 
pative term from which the 2D Euler equation is 


recovered for v = 0. The physics and solutions of the 
2D Euler and Navier-Stokes equations have been stud¬ 
ied extensively theoretically, numerically and in the lab¬ 
oratory, in the frame work of 2D hydrodynamic turbu¬ 
lence (see reviews bv iKraichnan fc MontgomervI 1 19801: 
iTabelin j I2QQ2I: iBoffetta fc Eckd I2Q12D . Unlike the 3D 
case, it has been shown that given a smooth initial condi¬ 
tion uo(x, y) at time t = 0 the 2D Euler equation admits 
a unique and regular solution at t > 0, i.e. , no finite time 
singularity develops (iRose fc Snleml 119781: iCheminI 119931: 


iBertozzi fc Constantinlll994 iMaida fc Bertozzill2QQl[). 

In 2D in addition to energy^ also mean-square vortic¬ 
ity (enstrophy) is conserved. The coupled conservation 
constraints have a strong impact on the dynamics that 
differs considerably from its 3D hydrodynamic counter¬ 
part and the corresponding magnetohydro dynamic cases. 
In particular, indicating the Energy with E = (1/2) (u^) 
(the integrated square velocity), the enstrophy with U = 
(1/2)and the palinstrophy with P = (1/2)(| Vcjp), 
the following energy and enstrophy conservation equa¬ 
tions are obtained from the 2D Navier-Stokes equations 
(e.g., IBoffetta fc Eckel[2Q12[ ): 


dt 


= —2uQ = —ey{t)^ 


dn 

dt 


= -2nP. 


( 11 ) 


Since all quantities (E, U, P, and u) are positively de¬ 
fined, it follows that U can at most decrease. Therefore 
the energy dissipation rate vanishes as viscosity tends 
to zero: 

( 12 ) 


lim Cy = 0. 


This result strongly differs from the 3D case where, 
i n the K41 phenomenology introduced by iKolmogorovI 
(ITmI . for a sufficiently small viscosity the energy dissi¬ 
pation rate is approximately constant Cy ^ const and 
independent from viscosity. This dissipa t ive an omaly 
in 3D was first pointed out by iTavlorl (1193511, and 
later confirmed in la boratory exper i ments (|Drvdenlll943l: 
iSreenivasanI 11984 iPearson et al.l 1200211 and by hy¬ 
drodynamic numer ical simulations (|Sreenivasa aiiM 
iKaneda et al.ll200^ . 

Thus, in contrast to the 3D case. Equation ([12]) implies 
that for small viscosities 2D turbulence is essentially un¬ 
able to dissipate energy at small scales. The viscous sink 
of energy is missing, because at any given time during 
the decay of an initial large-scale velocity field, for a suf¬ 
ficiently small value of z/, the dissipation is arbitrarily 
small. 

Therefore in two dimensions t here cannot be a d irect 
energy cascade. As proposed bv iKraichnanI (|1967[ ). en¬ 
ergy must flow toward the larger scales through an in¬ 
verse easeade. Thus during their evolution vortices (ve¬ 
locity eddies) acquire increasingly larger scale s, while it 
is en strophy that develops a direet cascade (|Batchelorl 
|1969[ ) with vorticity acquiring smaller scales. As can be 
seen from Equation the convective derivative of uj 
vanishes as z/ tends to zero, i.e., vorticity is constant in 
time at a point that moves with the fluid. The direct 
enstrophy cascade implies that an initial patch of vortic¬ 
ity gets stretched in time to form filamentary structures, 
so that while uj is convected with the fluid its gradient 
increases. 

Most numerical and laboratory experiments include a 
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small viscosity, but the regularity of the solutions of the 
2D Euler equation and the vanishing of the energy dissi¬ 
pation rate as viscosity tends to zero for the 2D Navier- 
Stokes equation allow to establish a clear connection be¬ 
tween the solutions of the 2D Euler and Navier-Stokes 
equations. Namely, the time evolution of the solutions 
of the same decay problem for the 2D Navier-Stokes and 
Euler equations is the same until dissipation sets in, i.e., 
until sufficiently small scales are formed and dissipation 
occurs in the Navier-Stokes case. 

The dual cascade picture has been investigated and 
supported by 2D Navier-Stokes simu l ations of forced 
turbulence (e.g., iHossain et al.l 119831: lErisch fc Sulem 
1984 : iSommeriS 119861: iTabelinj I2QQ2I: iBoffetta fc Ecke 
2Q12f ). of the decay of an initial condition consist¬ 
ing o f large-scale vorti c es (iMatthaeus fc MontgomervI 
I198QI: IMatthaeus et al.l Il991f ). and by laboratory 
experin ients of 2D decays performed with soap 
films. (jGharib fc Derangol 11989 : iBelmonte et al.l 119991: 
iGreffier et al.ll2QQ2l: iRivera et al. [20031). 

The decaying (initial value problem, or ‘relaxation’) 
case is pa rticu larly rele vant to the analysis carrie d out 
in Section l3Tl Recently iMininni fc Pouqu^ (j2Q13[ ) have 
performed a number of numerical simulations of the de¬ 
cay of an initial condition with energy in a narrow band 
of wavenumbers (with length-scale and velocity ui)^ 
with sufficient resolution to allow the development of 
both the inverse energy and direct enstrophy cascades. 
Since a decay is inherently non-steady, in order to com¬ 
pare with Kolmogorov K41 phenomenology, they per¬ 
form several simulations where the initial condition has 
the same velocity rms but different random ampli¬ 
tudes. In this way different realizations are obtained, al¬ 
lowing them to perform ensamble averages that smooth 
out the fluctuations in a single realization, and can then 
be compared more straightforwardly with the original 
K41 phenomenology (that as a matter of fact uses en¬ 
samble averages). Indeed a elear dual easeade is identified 
also in the ease of a deeay. In particular at wavenumbers 
smaller than the wavenumber of the initial condition, the 
peak of energy moves toward smaller wavenumbers with 
time, and energy develops an Ek ^ spectrum, fol¬ 

lowing K41 phenomenology (that does not depend on 
the direction of the cas cade - inverse or direct - e.g., 
see lR.ose fc SulemI 119781 ). The characteristic dynamic 
timescale is given by the eddy turnover time 

te = —, (13) 

U£ 

where i and ui initially are the length-scale and velocity 
of the initial condition. 

In K41 phenomenology the eddy turnover time m is 
the typical timescale for an eddy of size ^ ^ to undergo 
a significant distortion due to the relative motion of its 
components, thus transferring (in the 2D case) its en¬ 
ergy at larger scales, as schematically shown in Eigure[T] 
(left panel). The dimensional analysis of Equation ^ 
shows that is also the timescale over which an initial 
condition uq with energy at scales ^ £ and = ui 

undergoes a significant distortion with ^ ui. 

Eurthermore scales are defined as logarithmic bands of 
wavenumbers, e.g., kn G (2’^,2’^+^], with n G N and 
scale in ^ i± I kn = i± (the index n will be dropped 


hereafter). I n fact a single waven umber cannot repre¬ 
sent a scale (jAluie fc EvinkI l2QlQf ) since, from the un¬ 
certainty principle, the associated Eourier mode is de¬ 
localized in space and does not give rise to a localized 
physical structure such as an eddy^ the building block of 
K41 phenomenology. Consequently in the time interval 
ti [Equation i lT7|) / the energy of the system is transferred 
approximately from the seale i to the larger seale of dou¬ 
ble size 2£. 

Due to the structure of the absolut e statistical equi¬ 
libria of the ideal truncated sy stem (jKraichnanI 119671: 
iKraichnan fc Mont gomer"^ 119801 ) the development of an 
inverse cascade is expected in both the 2D dissipative 
(Navier-Stokes) and ideal (Euler) cases. Generally the 
time evolution of dissipative and ideal systems displays 
overall similar dynamics until energy reaches the small 
scales, when respectively thermalization and dissipation 
sets in. This is observed even when a direct energy cas¬ 
cade occurs, such as in 2D and 3D MHD, although spec¬ 
tral indices are observed to deviate from K41 in t he ideal 
cases (e.g., IWan et ^120091: iBrachet et al.ll2Q13[ ). Since 
the dynamics of the 2D Euler system may depart from 
K41 phenomenology, that has been developed and stud¬ 
ied for the dissipative case, in the following the eddy 
turnover time ti is then used as an estimate for the dy¬ 
namical timescale of the 2D Euler equation solutions. 
This is also justified by the dimensional analysis of an ini¬ 
tial condition consisting of vortices of scale i and velocity 
U£^ that indicates t£ ^ ilu£ as the order of magnitude of 
the dynamically relevant timescale. 

3.1. Strueture of Redueed MHD Equilibria and 
Dynamies 

This phenomenology can be applied to the reduced 
MHD equilibria (Equation [8]) using the mapping (fTQ|) . 
as schematically shown in Eigure [TJ Eddies^ that in 
the hydrodynamic case are velocity vortices, correspond 
now to magnetie islands. Then, given a magnetic field 
b6d(x, y) at the boundary z = 0 with energy at scales ^ i 
(i.e., a field structured in magnetic islands of scale ^ i)^ 
the unique and regular equilibrium solution heq{x^y^z) 
with heq{x^y^z = 0) = h^ix^y) is characterized by an 
increasingly stronger inverse easeade in the x-y planes 
for higher values of 2 : (see Eigure [H right panel), with 
the orthogonal magnetic field length-scale £ getting pro¬ 
gressively larger up to doubling its value in the plane 

(14) 

Obd 

the analog of the eddy turnover time, derived from Equa¬ 
tion m using the mapping m- If the magnetic field 
is characterized by a scale of order ^ at 2 : = 0 it will have 
its energy at scale ^ 2^ at 2 : = 2 :^, corresponding respec¬ 
tively to magnetic islands of scales £ and 2£ in the x-y 
planes 2 : = 0 and z = Z£. 

Therefore the equilibrium solution heq{x^y^ z) is gen¬ 
erally asymmetrie in the axial direetion z, but it ean be 
almost invariant or have strong variations, depending on 
the relative value of the length-seale Z£ eompared to the 
loop length Lz. As long as the variation scale is larger 
than the loop length {z£ > Lz) the field has a weak vari¬ 
ation along 2 ), but for increasingly smaller values of Z£ 
(z£ < Lz) the field becomes progressively more asymmet¬ 
ric along z due to the inverse cascade. Eixed the value 
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2D Euler equation 


Reduced MHD equilibria 



Figure 1. Schematic of the inverse energy cascade for the solutions of the 2D Euler equation (left panel) and for reduced MHD equilibria 
(right panel). In the hydrodynamic case (left panel) an eddy of size i doubles its size in the eddy-turnover time t£ ^ ijui (Equation IE]). 
For the structure of reduced MHD equilibria (right panel) this implies (via the mapping t ^ z, u —)■ b/Bo given by Equation (|10I) I that an 
equilibrium solution with magnetic islands of transverse scale i in the plane z = 0 will acquire an axial asymmetry along the ^-direction 
characterized by the variation length-scale Z£ ~ Boijbi^d (Equation [H]), since magnetic islands in the plane z — zg have an approximately 
double transverse scale (~ 2^). 


of the guide field Bq and the scale of the magnetic field 

the axial variation scale zg is inversely proportional to 
the magnetic field intensity at the boundary Bm- quasi- 
invariant equilibria are then obtained for weak magnetic 
fields while increasingly stronger fields result in more 
asymmetric equilibria. 

A strong asymmetry of the magnetic field b along z 
is generally not eompatiBle with the dynamical solutions 
of the reduced MHD equations ([I])-(|2j), as confirmed by 
nonlinear simulations (for a more detailed discussion of 
the topics summarized in the present and next paragraph 
see iRappazzo et al.l l2QQ8[ ) . The derivatives along z ap¬ 
pear only in linear terms in Eqs. o-®. Introducing the 
Elsasser variables = u ± b, and neglecting nonlinear 
and diffusive terms, the remaining linear terms yield the 
two wave equations: 

with V • = 0. (15) 

Thus fluctuations propagate along z at the Alfven speed 
Bq, and a strongly asymmetric field cannot be generated, 
particularly for the problem considered here with line- 
tying boundary conditions. Eor instance a velocity Uph 
at the boundary z = 0 implies the “reflection” condition 
z~ = —z+ -b 2uph (z“ propagates inward, z+ outward), 
i.e., generalized Alfven waves are injected and propagate 
in the computational box, and only the m = 0 mode does 
not propagate in the axial direction. 

Considering an initial condition with only the guide 
field B = BqBz^ vanishing orthogonal component b = 0, 
and a constant velocity Uph at the photospheric bound¬ 
ary that shuffles the magnetic field line footpoints, an or¬ 
thogonal component of the magnetic field invariant along 
2 : is generated and it grows linearly in time as 

h{x,y,z,t) = Uph{x,y) — , (16) 


where ta = LzIBq is the Alfven crossing time in the ax¬ 
ial direction. Strictly speaking higher modes are present 
depending on how the velocity Uph is turned on, but 
they represent a small contribution compared to Equa¬ 
tion m that is the strongly dominant term. 

The solution (Hbl) is obtained from Equations ([T^ (see 
IRappazzo et al.l I 2 QQ 8 I ). thus negleeting nonlinear terms 
in the reduced MHD e quations. This is just i fied a s long 
as B is small. In fact IRappazzo fc Parl^ (j2Q13[ ) have 
shown that for initial configurations with only the m = 0 
mode in the axial direction and a non-vanishing 2D or¬ 
thogonal Lorentz force component (with corresponding 
term b • Vj 7 ^ 0 in Equation [ 8 ]), nonlinearity is strongly 
suppressed, i.e., the nonlinear terms ean Be negleeted, 
for magnetie fields with intensities Below the magnetie 
threshold B ^ iBojLz. The magnetic field decays only 
for larger magnetic fields B > iBojLz, while for smaller 
values the system does not form significant current sheets 
and energy is not dissipated. 

The intensity threshold B ^ IBq/Lz corresponds to a 
eritieal equilibrium variation length-seale of the order of 
the loop length zg ^ Lz (Equation [H]). Since fields with 
only m = 0 are invariant along ^ the correspondent equi¬ 
libria length-scale zg can be computed with Bm = B. Eor 
B < IBojLz the corresponding equilibrium, i.e., the equi¬ 
librium solution computed with Bm = 5 at the boundary, 
has Zg > Lz and it is quasi-invariant along z. Conse¬ 
quently magnetic fields with m = 0 and intensity below 
the intensity threshold B ^ £Bq /Lz are very close to their 
corresponding equilibrium solution. It is sueh elose prox¬ 
imity to an equiliBrium that suppresses nonlinearity, that 
at equilibrium is indeed entirely depleted. 

The emerging phenomenology for the dynamics of ini¬ 
tially straight axial field lines shuffled by a velocity 
field Uph constant or slowly changing in time (so that 
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the induced magnetic field is quasi-invariant along z) is 
therefore the following: since at first the induced mag¬ 
netic field is small, the associated variation scale is large 
Lz, nonlinearities are suppressed and the magnetic 
field grows as in Equation until the variation length- 
scale becomes smaller than the loop length Z£ < Lz, when 
nonlinearity can develop leading to the formation of cur¬ 
rent sheets. 

The decay of initial configurations with a parallel m=0 
mode are relevant for slow photospheric motions. But in 
general the system has three characteristic timescales: 

1. the surface convective timescale Tsc ^ 
essentially the typical lifetime of a granule, approx¬ 
imately given by the ratio of half its length-scale I sc 
over the photospheric velocity Uph (typical values 
for the Sun are I sc ^ 10^ km, Uph ^ Ikm/s, and 
Tsc ^ 5 - 8m), 

2. the Alfven crossing time ta = LzjBo, where Lz 
is the loop axial length and Bq the Alfven length 
associated to the guide field, and 

3. the nonlinear timescale Tnr that is investigated 
theoretically and numerically. 

For typical X-ray bright loops Lz ^ 4:0 x 10 ^ km and 
Bq 2 X lO^km/s, therefore the Alfven crossing time 
Ta = LzIBq 20 s is much smaller than the photo¬ 
spheric timescale, with taItsc ^ 0.04. For these loops 
photospheric motions are then characterized by a low 
frequency, i.e., 

Tsc > ta, ( 17 ) 

and a constant {zero frequency, Tsc = oc) photospheric 
velocity can be a good approximation, since such slow 
variation of photospheric motions {tsd^A ^ 25) intro¬ 
duces only wavelenghts much longer than the loop length 
itself along z, and the resulting magnetic field (pT|) can 
be considered invariant along z. 

Nevertheless this condition can break down for longer 
solar coronal loops and for loops on other active stars 
with outer convective envelopes and magnetized coro- 
nae, that exhibit broad variations in magneti c field 
intensity and topologies (jPonati fc Landstreetl I2QQ9I: 
iReinersI l2Q12l), and photospheric m otion properties 
(|Ludwig et al.l l2QQ2l : iBeeck et al.l l2Q13[ ) , for which Tsc ^ 
ta^ or Tsc < ta- In this case the resulting magnetic field 
will have a more complex expression than Equation ([T6]) , 
and higher modes along 2 : will be present and contribute 
increasingly more, the faster the convective timescale Tsc 
compared to the Alfven crossing time ta- 

Therefore the structure of the magnetic field induced in 
coronal loops by photospheric granulation will be domi¬ 
nated by the m = 0 mode along 2 : (i.e., the field is quasi¬ 
invariant along z) for the typical X-ray bright solar loops 
for which Tsc » ta, while for longer loops (including 
loops on other active stars) higher modes (m > 1) will 
be increasingly more important the smaller the timescale 
ratio TscIta < 1- 

In both cases the variation scale Z£ ^ IBofbM (Equa¬ 
tion [14]) measures the asymmetry of the equilibrium so¬ 
lution. In both cases, even in presence of higher modes 
m > 1, the dynamical solutions of the reduced MHD 
equations will not be asymmetric along 2 :, in contrast to 
the equilibria with Z£ < Lz. 


Table 1 

Simulations Summary 


Run 

z-modes 

Bo 

numerical grid 

Re 4 

A 

m=0 

10® 

2D: 5122 

1019 

B 

m=0 

10® 

3D: 5122 252 

1019 

C 

m=l 

10^ 

3D: 5122 252 

1019 

D 

m=0-4 

10® 

3D: 5122 252 

IOI 9 


Note. — The second column indicates the parallel Fourier 
modes used in the initial magnetic field (Equation [3]). So is 
the axial Alfven velocity (same for all runs). The numerical 
grid Ux X Uy X Uz is indicated in the fourth column, it is three- 
dimensional for all runs except for run A that is two-dimensional 
(with grid Ux X Uy). The last column indicates the value of the 
hyperdiffusion coefficient Re^ = —l/i >'4 = — 1/774 used in Equa¬ 
tions 0 - 0 . As described in the text, each run B-D is a collec¬ 
tive label for a set of simulations carried out with the parameters 
indicated in the table, but with initial conditions (Equation [3]) dif¬ 
fering for the values of the ratio 60 /So, the rms of the orthogonal 
magnetic field over the guide field intensity. 


In the following sections the dynamics of configurations 
with different initial conditions, with a single mode along 
z {m = 0 and m = I), and with all modes 0 < m < 4 
excited, are investigated through numerical simulations 
(see Tabled]). 

The 2D photospheric velocity Uph models photospheric 
motions. These have large scales (with £ r\j 10 ^ km) 
and are disordered, since they originate from turbulent 
convection. Thus generally the magnetic field b will 
have a non-vanishing 2D Lorentz force component, i.e., 
b • Vj 7 ^ 0 , since the opposite would imply that j should 
be constant on the field lines of b, a condition too sym¬ 
metric to apply to turbulent convection (this could be re¬ 
alized, e.g., with an exactly ID magnetic shear along one 
direction, or with perfectly circular field lines). There¬ 
fore in all simulations presented here the 2D Lorentz force 
component of the magnetic field does not vanish, i.e., 
b • Vj 7 ^ 0 . 

Notice that when b • Vj = 0 we obtain dzj = 0 from 
the equilibrium Equation ^ , hence in this case the vari¬ 
ation length-scale is formally infinite Z£ = 00 and it is 
always larger than the loop length. The inverse cascade 
picture of the 2D Euler equation described in section [3] 
is indeed valid only for initial conditions that are out of 
equilibrium. If u • Vcj = 0 in Equation (|9]) then there is 
no time evolution and formally T£ = oc (Equation [l3]). 
Reduced MHD equilibria with dzj = 0 (or equivalently 
b • Vj = 0 ) are those apt to describe classic linear insta¬ 
bilities (such as kink instabilities), and their dynamical 
accessibility will be discussed in the following sections. 

4. RESULTS 

The results of the numerical simulations are presented 
in this section. All simulations, with the exception of 
Run A, implement line-tying boundary conditions with 
field line footpoints rooted in a motionless plasma in the 
photospheric-mimicking planes 2 : = 0 and z = Lz. Run A 
implements periodic boundary conditions in the axial di¬ 
rection 2 :, and since its initial condition is invariant along 
2 : the dynamics will not introduce any variation along 
along this direction {dz = 0 ) and the simulation is re- 
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stricted to a 2D plane. In the orthogonal directions (x-y) 
all runs use periodic boundary conditions. 

All simulations consider the decay (or equivalently re¬ 
laxation) of an initial magnetic configuration ma de of 
large-scale Fourier modes (as described in Section \2Jh 
with non-vanishing 2D orthogonal Lorentz force compo¬ 
nent (b • Vj 7 ^ 0). The guide field intensity is Bq = 1000 
for all simulations, corresponding to an Alfven velocity 
of 1000 km/s, a typical value for solar coronal loops. 

Dissipative simulations with initial conditio ns m ade of 
sing le pa rallel modes are describ ed in Sections l4T] (m=0) 
and 14.21 fm=lb while in Section lT3l the initial condition 
has all modes with 0 < m < 4. A summary of the 
simulations is shown in Table [H 

The initial magnetic fields used in these decaying sim¬ 
ulations can be considered as “snapshots” of the coronal 
magnetic field at different stages of its evolution, partic¬ 
ularly for the problem of an initially straight field shuf- 
fied at its footpoints by photospheric motions [see Equa¬ 
tion ([T6|) ]. The different parallel modes are related to the 
photospheric motions frequency, as discussed in the fol¬ 
lowing sections. Additionally the relaxation of magnetic 
fields is a to pic of gener al and broad interest to solar 
physics fe.g.. |Pri^l2Q14f) . 

4.1. Runs A-B: single mode m=0 

This section presents an extended analysis of the dis¬ 
sipativ e simulations carried out by iRappazzo fc Parkerl 
(|2Q13f) that investigate the decay (or equivalently relax¬ 
ation) of initial magnetic configurations with an out-of- 
equilibrium orthogonal magnetic field and parallel mode 
m=0 (runs A-B). 

The initial magnetic field is invariant along ^ since only 
the mode m=0 is present. It has the same structure for 
all runs A-B, with same topology for the magnetic field 
lines of b, as shown in Figure |4] at time t = 0, but dif¬ 
ferent values of the orthogonal field intensity rms ho (the 
proportionality factor in equation [3]) are used. Respect 
to the guide field Bo (same for all runs) the intensity is 
ho = 0.1 Ro foi* run A, while run B comprises a set of 
simulations performed with same parameters except ho 
that spans the range 0.01 < ho/Bo < 0.1. 

Runs A and B differ for the boundary conditions along 
z, periodie for run A and line-tied for runs B. In the 
periodic case the invariance along ^ is preserved during 
the time evolution, therefore for run A Equations 0-® 
are integrated in a 2D x-y plane with 9;^ = 0. 

Periodicity along 2 ; is appropriate as a local approxi¬ 
mation for the central part of very long loops, where line- 
tying at the boundary has a weak inffuence. In partic¬ 
ular line-tied forced simulations, with a velocity field at 
the photospheric-mimicking boundaries, have shown that 
dynamics resembles those of periodic sin iulations for low 
values of t he ratio fc = iscBo/LzUph (jRappazzo et al.l 
I2QQ7L [ 200 ^ . Fixed the scale of granular cells I sc and the 
photospheric velocity rms Uph^ line-tying has a weaker 
impact on longer loops (with larger Lz) or loops with 
weaker guide fields (smaller Bo). 

Since the 2 D run A is started with an out-of- 
equilibrium magnetic field (b • Vj 7 ^ 0 ), it is therefore 
akin to 2D turbulence decay simulation s (jHossain et al.l 
I 1995 I: iGaltier et al.lll99^ lBiskampll2QQ3f ) that use similar 
initial conditions for the magnetic field, but addition¬ 
ally have an initial velocity field in equipartition. The 



0 10 20 30 t/t^ 40 50 60 

Figure 2. Runs A-B (m=0): Total energy vs. time for line-tied 
simulations with different values of 60/^0 (runs B) and the 2D sim¬ 
ulation (run A, with ho/Bo = 10%). The inset shows in logarithmic 
scale total energy normalized with its inital value. 


dynamics are nevertheless similar and in the 2D case to¬ 
tal energy decays approximately as E oc t~^ (Figure O 
inset), even if the initial velocity vanishes everywhere. 
Until time t ^ 0.4 total energy is conserved, but the 
non-vanishing Lorentz force transfers ^ 15% of magnetic 
energy Em into kinetic energy henceforth leading 
to the formation of small-scales and current sheets, dis¬ 
sipating ^ 90% of the initial magnetic energy, while ki¬ 
netic energy remains much smaller tha n magnetic energy 
throughout (|RaDDazzo fc Parkerll2013h . 

In Fourier space magnetic energy has initially only 
large-scale perpendicular modes k=3 and 4 (Figure [31 
top panel). As nonlinearity develops energy is progres¬ 
sively transferred toward both small (direet) and large 
scales {inverse easeade). In physical space the direct cas¬ 
cade gives rise to current sheets (Figure (H top row) that 
enable dissipation through magnetic reconnection. As 
dissipation peaks at t ^ 1.7 the spectrum extends 
fully toward high wave-numbers exhibiting an approxi¬ 
mate power-law. At the same time a substantial 

fraction of total energy has already been transferred to 
large-scale modes k=l and 2 through the inverse cascade, 
that in physical space corresponds to the eoaleseenee of 
magnetic islands as magnetic reconnection occurs. As 
dynamics proceeds the energy transferred at the small 
scales is dissipated leading to the disappearance of cur¬ 
rent sheets and small-scales, so that the system finally 
relaxes to a state with energy mostly at large scales 
(particularly in mode k=l), corresponding in physical 
space to large-scale magnetic islands and large-scale cur¬ 
rent layers (Figure 01 top row, t ^ 2 10 r^). The w hole 
process is akin to Taylor relaxation (|Tavlorl [l986D and 
self-organi zation leading to the formation of large-scale 
structures (|Hasegawalll985f ). 

In the 2 D case (run A) solutions differing only for the 
value of ho (the rms of the initial magnetic field in Equa¬ 
tion [3]) have a self-similar structure. Indicating with 
' 0 o(x, t) and (/ 9 o(x, t) the solution of Equations ([I])-© 
with initial magnetic field rms 60 , then the solutions with 
h'o = cr 6 o, and same random amplitudes in Equation ([3]), 















































9 


are given b£l: 

Ip' (x, t) = a ipo (x, at), (fi' (x, t) = a ipo (x, at) , (18) 

as can be verified by direct substitution. Consequently 
all these solutions have a similar structure and their time 
evolution differs only for the scaling factor a. In partic¬ 
ular if current sheets form for a certain value of 6o, they 
will always do for any value of 6o at scaled times. Anal¬ 
ogously energy will exhibit a power-law decay with the 
same exponent because E'(t) = a‘^Eo{at), implying that 
if E{t) oc then E'{t) oc 

When the same initial condition is used with line-tying 
boundary conditions, the system is no longer invariant 
along 2 ), as now the velocity must vanish at the top and 
bottom plates z = and z = L, therefore the velocity 
cannot develop uniformly along z as in the periodic case. 

The time evolution of total energy for line-tied simu¬ 
lations with different values of bo is shown in Figure [2l 
While the dynamics of the system with bo/Bo = 10% is 
similar to the 2D case with energy dissipating ^ 84% of 
its initial value, the behavior is increasingly different for 
lower values of 6o, with progressively less energy getting 
dissipated. For bo/Bo ^ 3% no significant energy dissi¬ 
pation nor decay are observed. Additionally, also for the 
decaying cases their dynamics are strongly suppressed 
once energy crosses this threshold. As shown in Fig. [2] no 
energy decay is observed below ^ 5 x 10^, correspond¬ 
ing to a ratio /Bo ^ 3%. The inset in Fig.[2]shows 

that energy decays with different power-law indices for 
lower values of bo/Bo, hence time self-similarity is lost 
and the impact of line-tying on the dynamics is more 
complex than a simple delay as in the 2D case (Equa¬ 
tions [H]). 

Magnetic energy spectra (integrated along 2 ), Figure [3|) 
show similar dynamics for the 2D {top panel) and the 
3D case with bo/Bo = 10% {middle panel). In par¬ 
ticular the spectra are fully extended toward the small 
scales as dissipation occurs, corresponding to the forma¬ 
tion and dissipation of current sheets in physical space 
(Figure m seeond row). While similar behavior is ob¬ 
served for magnetic field intensities bo > 3%, below this 
threshold the spectra do not extend to the high wave- 
numbers where energy gets dissipated (Figure [3l bottom 
panel, bo/Bo = 2%), i.e., eurrent sheets do not thin be¬ 
low the eritieal diffusive thiekness that allows magnetie 
reeonneetion and energy dissipation to oeeur. As shown 
in Figure m at the peak of dissipation {eentral eolumn) 
both current maxima and the number of current sheets 
decrease for smaller bo/Bo, and for bo/Bo = 2% no cur¬ 
rent sheets are formed, but only a few ripples are visible 
in the magnetic field (enhanced in the current) and are 
eventually dissipated on long timescales. 

Furthermore, spectra show that the inverse energy 
cascade decreases from the 2D to the 3D case with 
bo/Bo = 10%, in fact while in the asymptotic state of the 
2D run {t = 210 r^) most of the energy is in the k = 1 
mode and higher modes are much smaller, in the 3D case 
the mode k = 2 has the higher value. For 3D simulations 
with smaller bo/Bo the inverse cascade is progressively 

^ Strictly speaking these self-similar solutions would require the 
Reynolds number to scale as R' = aR, but in the high-Reynolds 
regime the solutio ns of decaying turbulence do not d epend on the 
Reynolds number (IBiskamDll2003l : iGaltier et al.lll997| b 





Figure 3. Runs A-B (m=0): Magnetic energy spectra (integrated 
along z) at selected times for the 2D run A with bo/Bo = 10% {top 
panel)^ and for run B simulations with ho/Bo = 10% (middle) and 
bo/Bo = 2% (bottom). Energy is normalized with its initial value 
at time t=0, when energy is present only at perpendicular modes 
k = 3 and 4 (diamond symbol). Spectra at the time of maximum 
dissipation for each simulation are drawn in a continuous line. 


fainter and does not occur for bo/Bo ^ 3%, as shown in 
Figure [3] {bottom panel) for bo/Bo = 2%, where no sig¬ 
nificant energy is found in perpendicular modes smaller 
{k < 2) than those present at time t = 0 {k = 3, 4). 

In physical space (Figure 2]) the inverse cascade cor¬ 
responds to larger-scale magnetic islands in the asymp¬ 
totic state. Since a larger fraction of magnetic energy 
is dissipated for higher values of bo/Bo and in the 2D 
case, more magnetic fiux is reconnected, thus leading to 
increased eoaleseenee and larger magnetic islands. Ad¬ 
ditionally the field line topology in the relaxed state is 
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Figure 4. Runs A-B (m=0): Magnetic field lines of the orthogonal magnetic field component b and current density j at selected times. 
Top row shows snapshots from the 2D simulation (run A), while snapshots from three 3D simulations (run B) with different ratios bo/Bo 
are shown in the remaining rows (in this case the mid-plane 2 ; = 5 is considered). The first column shows the initial condition at t = 0, 
same for all except for the ratio bo/Bo, the second column shows the fields at the time of maximum dissipation, while the third column 
shows the fields at a later time when the fields have relaxed and little if any energy dissipation occurs. 


substantially unaltered respect to the initial condition 
for bo/Bo < 3% (cf. the first and last columns in Fig¬ 
ure 12 ), with higher variations for higher values of bo/Bo 
and most of all in the 2D case. 

A key difference distinguishes the 2D and 3D asymp¬ 
totic topologies. Although all simulations are started 
with a non-vanishing 2D Lorentz force component (b • 
7 ^ 0)5 fhe 2D simulation relaxes to an approximate 
equilibrium with b • Vj = 0, but all 3D simulations 


relax to an orthogonal field with b • Vj 7 ^ 0 , regard¬ 
less of how much energy is dissipated during the decay 
(none for bo/Bo = 2%, and an 84% energy decay for the 
bo/Bo = 10% run B). Further analysis is presented in the 
following to understand the nature of these equilibria and 
how they are approached. 

In the 2D case the equilibrium condition (| 8 |) becomes 
simply b • Vj = 0. This requires the current density j to 
be constant along the field lines of b (or equivalently that 
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the isosurfaces of j are also isosurfaces of the magnetic 
potential t/^). This condition is generally satisfied only in 
highly symmetric configurations such as one-dimensional 
magnetic shears, e.g., b = /(^) or rotationally invari¬ 
ant fields as b = /(r) bq (/ is a generic function in carte¬ 
sian or cylindrical coordinates). In the 2D case the field 
lines are mostly circular in the asymptotic state (Fig¬ 
ure m and the isosurfaces of the current density and of 
the magnetic potential (the field lines of b) overlap. In 
the 3D case the full equilibrium equation (jHj) has to be 
considered, and the fact that for run B with bo/Bo < 3% 
no significant dynamics occurs with the initial orthogo¬ 
nal magnetic field essentially unaltered even though its 
2D Lorentz force does not vanish, i.e., b- Vj 7 ^ 0, implies 
that dzj increases only slightly from its initial vanishing 
value. 

The dynamical approach of the system to equilibrium 
is further investigated with the probability density func¬ 
tions (PDFs) of the equilibrium equation terms. These 
PDFs are histograms of the quantities of interest nor¬ 
malized so that the resulting function f{q) multiplied by 
the bin size Aq gives the fraction of points in the compu¬ 
tational box where the specific quantity q has its value 
in the interval [q — Ag/ 2 ,g + Ag/ 2 ], and the integral 
f ^Qf{Q) = 1- PDFs have been computed for both the 
left and right hand terms in the equilibrium equation (j8j) , 
and for their difference, indicated with 


This is a three-dimensional scalar function that measures 
how far the system is from equilibrium locally at point x, 
vanishing at equilibrium and with higher absolute values 
the larger the departure from equilibrium. The PDFs 
have been computed for all the simulations with m = 0 
(for the 2D case only the term b • Vj / Bq is present), and 
they all exhibit similar behavior. All these PDFs have 
vanishing averages, therefore their standard deviations 
can be defined as: 


= 




^eq 


9 \ 1/2 

dj , b-Vjy\ 
dz^ Bo ) / ’ 


( 21 ) 


respectively for the first and second terms and the whole 
Equation m, labeled as parallel {cz), orthogonal (a^) 
and total {cr^. 

In Figure [5] {top panel) the PDFs of the equilibrium 
function Eq{yi) are shown in a semi-log plot at selected 
times for run B with Bq/Bq = 10%. The abscissa is 
rescaled with the standard deviation of the PDFs to im¬ 
prove visualization, since they exhibit large variations in 
time. The PDF of Eq is generally super-Gaussian (with 
a peak around zero and “tails” farther out), particularly 
close to maximum dissipation time {t = 2ta)^ however 
the central part appears closer to a Gaussian distribution 
at later times {t = lOr^), and particularly in the final 
asymptotic stage {t = 1018.2 r^). Although long tails 
are present at most times, in all cases ^ 95% of points 
lies within two standard deviations from zero (from 94% 




Figure 5. Runs A-B (m=0): Probability density functions 
(PDFs) of the reduced MHD equilibrium equation ^ at selected 
times for the run B with bo/Bo = 10% (top panel), shown in a 
semi-log plot. To accommodate the large variation of the standard 
deviation the abscissa displays the values normalized with the stan¬ 
dard deviation at the corresponding time. The remaining panels 
display in logarithmic scale the standard deviations as a function 
of time for the 2D simulation (bottom panel) and the 3D runs B 
with bo/Bo = 2% and 10% (middle panels). For the 3D cases the 
orthogonal (cr^), parallel (az) and total (cTeq) standard deviations 
are shown. 
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at t = 0 to 96% at t = 10 r^). 

An appropriate quantitative measure of the distance of 
the system from equilibrium is given by the standard de¬ 
viation of Eq^ shown in the mid panels of Figure [5] along 
with (Tz and for runs B with bo/Bo = 2% and 10%. 
In this case, since their averages vanish, the standard de¬ 
viations are also the rms of the considered quantities. At 
time t = 0 the derivative along 2 : of j vanishes {dzj = 0) 
so that initially az = 0 for both runs, while = 25.32 
and 633.09 respectively. But already after one Alfven 
time ta they both increase (substantially only for the 
run with bo/Bo = 10%) reaching similar values az ^ 
and subsequently continue to be very close while their 
values decrease. The rms of the equilibrium function Eq^ 
the standard deviation deg, decreases with time, and it 
is asymptotically smaller than and a_\_. As shown in 
Figure 0 all standard deviations decrease asymptotically 
like a power-law with az ^ cr± (X where spectral 
indices are respectively a = 0.17 and 0.038 for the runs 
with bo/Bo = 10% and 2%, while the rms of Eq decays 
at a faster rate with aeq oc where /3 = 1.45 and 1.25. 

This implies that while the rms of dzj and b • Vj /Bo 
have about the same value and remain approximately 
constant in the asymptotic state (when their power-law 
decay occurs), equilibrium is approached as the two terms 
in Equation ([19]) balance each other progressively more 
throughout the computational box, thus leading to the 
rapid decrease of aeq. 

The initial increase of the standard deviations is larger 
for bo/Bo = 10% than for the 2% case, since for bo/Bo = 
10% the system is farther out of equilibrium in the begin¬ 
ning, with strong currents forming quickly and maximum 
dissipation rate occurring at t ^ 1.7 when also stan¬ 
dard deviations approximately peak. Their subsequent 
decrease up to t ^ 20 is enhanced by the strong de¬ 
cay of magnetic energy and progressive disappearance of 
current sheets, before approaching the asymptotic stage 
around t ^ 200 Although the standard deviations 
for run B with bo/Bo = 2% have similar behavior to 
the case with bo/Bo = 10%, their variations are much 
smaller, since the field just undergoes a slight adjust¬ 
ment, without any significant energy decay (Figure [2|), 
current sheets formation or significant change in the mag¬ 
netic field topology (see bottom row in Figure |4|). 

For the 2D run A there is no parallel standard devia¬ 
tion az. The rms of Eq, aeq^ is equal to a±, and simi¬ 
larly to runs B its time evolution follows the power-law 
aeq OC t~^ with p = 1.33, as shown in Figure [S] (bottom 
panel). It is worth mentioning that aeq initially grows 
larger respect to its initial value (= 633) for run B re¬ 
spect to run A with same bo/Bo = 10%. The reaction of 
the line-tied field sets the system further out of equilib¬ 
rium, an enhancement of nonlinearity that might favor 
the development of singularities in the line-tied system 
respect to the periodic case. 

Further insight into the dynamics is gained analyzing 
the probability density function of the cosine of the an¬ 
gle 0}) between b and Vj 


cos Ob = 


bVj 

|b||Vi|’ 


( 22 ) 


shown in Figure [6l At time t=0 the “2D perpendicu¬ 
lar” Lorentz force term does not vanish for all simula¬ 







Figure 6. Runs A-B (m=0): Probability density functions 
(PDFs) of the angle 6^ between the orthogonal magnetic field b 
and the current density gradient Vj at selected times for the 2 D 
run A with 60/^0 = 10% {top panel), and the 3D runs B with 
respectively bo/Bo = 10 % (middle) and bo/Bo = 2 % (bottom). 


tions (b • Vj % 0). Since the orthogonal component of 
the magnetic field differs for runs A-B only for a pro¬ 
portionality factor (Equation [3]), the PDE is the same 
for all runs A-B (e.g., see mid-panel in Eigure[6|). Al¬ 
though it is peaked around zero (corresponding to the 
2D equilibrium condition b • Vj = 0) it spreads out with 
significant values up to | cos^^j <1/2, corresponding to 
an approximate 60° angle around Ob = 90°. Eor the 2D 
run A with bo/Bo = 10% (top panel) the PDE initially 
spreads out during the strongest part of the nonlinear 
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stage (1 ^ t/rA ^ 10) during which ^ 90% of the initial 
energy is dissipated, but afterward peaks progressively 
more strongly around zero, corresponding to the equilib¬ 
rium condition for the orthogonal field b • Vj = 0, with 
respectively 71% and 96% of the grid points in the vol¬ 
ume in the region | cos 6^51 < 0.1, spanning an angle of 
^ 11.5° around 6^5 = 90°, at times t=29.1 and 911.6 ta- 
This confirms that the system approaches asymptotically 
an equilibrium with b • Vj = 0 corresponding in physi¬ 
cal space to increasingly circular field lines progressively 
more coincident with the isosurfaces of j, as shown in 
Figure m (top row). 

In contrast the picture is radically different for the line- 
tied simulations. As shown in the middle panel of Fig¬ 
ure [6l for run B with bo/Bo = 10% the PDF spreads 
increasingly further out during the time evolution^ flat¬ 
tening considerably already after one Alfven time and 
remaining flat throughout the subsequent evolution^ when 
^ 80% of magnetic energy is dissipated, and in the fol¬ 
lowing asymptotic regime, with peaks forming in corre¬ 
spondence of alignment between the two fields (^5 ^ 0°, 
180°). Therefore, in contrast with the periodic case (2D 
run A) the orthogonal magnetic field does not approach 
the asymptotic equilibrium with b • Vj = 0 (that would 
imply also dzj = 0 in the 3D case. Equation [8]), in¬ 
stead in the 3D line-tied case the orthogonal component 
of the magnetic field remains with a non-vanishing 
perpendiculaP^ Lorentz force term (b • Vj % 0). 

Furthermore, if the initial mag netic field intensity 
i s bel ow the threshold set out by iRappazzo fc Parl^ 
dMl), as shown in the bottom panel of Figure [6] for 
run B with bo/Bo = 2%, the PDF starts flattening out 
to some extent but then bounces back very close to its 
initial profile, corresponding to a slight readjustment of 
the magnetic field as shown in Figure H] (bottom row), 
with the orthogonal magnetic field preserving its non¬ 
vanishing perpendicular Lorentz force. 

The results of the numerical simulations analyzed in 
this section are consistent with the heuristi c phenomenol¬ 
ogy laid out by IRappazzo fc Parl^ (j2Q13f ) and the more 
refined analysis of the structures of the equilibria ex¬ 
pounded in Section [H The asymmetry along z of the 
solutions of the reduced MHD equilibrium equation ([8]) 
can be estimated with the axial variation length-scale 
Z£ iBo/b (Equation [14], see also EigureHj), where i is 
the perpendicular characteristic scale (in the x-y plane) 
of the magnetic field component b. 

As discussed in Section EH the dynamical solutions 
of the reduced MHD equations generally cannot 

exhibit strong asymmetries along 2;, in particular when 
driven from the photospheric boundaries. On the other 
hand the equilibria can be strongly asymmetric or quasi¬ 
invariant along z, depending on the relative value of the 
axial variation scale Z£ respect to the loop axial length 
Lz, with the critical length given by Z£ ^ Lz. Eor small 
values of b the axial variation scale Z£ is longer than the 
loop length Lz and the corresponding equilibrium solu¬ 
tion is quasi-invariant along while for larger values of 
b the axial scale Z£ is smaller than the loop length and 
the equilibrium is more asymmetric along 2; the larger 
the magnetic field intensity b. 

As shown in Eigure [2] no substantial energy is dissi¬ 
pated for 3D runs with bo/Bo ^ 3%, and just mini¬ 
mal dynamics occurs as the field slightly readjusts (Eig- 


ures[3]-[6|). Since the initial magnetic field (Equation [3]) 
has a perpendicular scale i ^ L±/3.S7 ^ 1/3.87 (the 
averaged wavenumber of the initial condition is 3.87 
and L± = 1), this threshold corresponds to a varia¬ 
tion length-scale for the initial magnetic field of about 
Z£ = iBo/bo ^ 100/(3x3.87), i.e., Z£ > Lz since Lz = 10. 

Therefore for bo/Bo ^ 3% the corresponding equilibria, 
computed from Equation [8] with heq{z=0) = bo(2:=0), 
as described in Section [H have a large variation length- 
scale Z£ > Lz and are therefore quasi-invariant along 
2;. Since the initial condition has only the parallel m=0 
mode it is invariant along 2;, and therefore very close 
to the corresponding equilibrium, so that nonlinearity is 
strongly depleted and only a slight readjustment of the 
field occurs, with no significant energy dissipation. On 
the contrary for initial conditions with bo/Bo > 3% the 
corresponding equilibria have smaller variation length- 
scales Z£ < Lz^ hence the equilibrium solution is strongly 
asymmetric along differing substantially from the ini¬ 
tial condition that is then necessarily out-of-equilibrium, 
as shown by the subsequent dynamics and energy dissi¬ 
pation. 

Eurthermore, as shown in Eigure [S] also in the cases 
when decay occurs, i.e., for initial conditions with 
bo/Bo > 3%, the magnetic field relaxes to a new equilib¬ 
rium that approximately satisfies the condition b/Bo ^ 
3% and Z£ ^ Lz. But while the asymptotic energy of 
the runs with bo/Bo = 4%, 5% and 6% are approxi¬ 
mately the same, corresponding to a ratio b/Bo ^ 3.3%, 
the run with bo/Bo = 10% relaxes to a slightly higher 
energy with b/Bo ^ 4%. On the other hand the run 
with bo/Bo = 10% has a stronger inverse cascade (Eig¬ 
ure 0), with significantly larger magnetic islands in the 
asymptotic regime (Eigure |4|) and average wavenumber 
^2.7 thus obtaining again Z£ ^ Lz. When a strong in¬ 
verse cascade occurs the formation of larger perpendic¬ 
ular scales {t) increases the value of the axial variation 
length-scale Z£ ^ iBo/b^ thus attaining the equilibrium 
condition Z£ ^ Lz with a larger value of 6, and conse¬ 
quently a smaller dissipation of energy. 

The critical variation length-scale Z£ origins from a 
balance of forces because the reduced MHD equilib¬ 
rium equation (j8|) represents a balance bet ween two force 
term s. In the reduced MHD limit (e.g., iMontgomervI 
Il982f ) the Lorentz force splits into two terms with com¬ 
ponents only in the orthogonal x-y planes: the “perpen¬ 
dicular” (b • Vb) and “parallel” {Bodzh) field line ten¬ 
sions (plus the pressure term, determined through the 
incompressibility condition). The first term represents 
the field line tension of the orthogonal component b, the 
only one present in the 2D limit. The second is an addi¬ 
tional tension term due to the presence of the guide field 
Ho, linked to the tension of the field lines of the total 
magnetic field Hoe;^ + b. An equilibrium is attained only 
when these two counteracting components of the Lorentz 
force ba lance each other satisfying Equation (j8|) . As out¬ 
lined in IRappazzo fc Parkerl (f2Q13f ) these two forces are 
of the same order of magnitude for the critical intensity 

(23) 

corresponding to the critical axial variati on l ength-scale 
Z£ ^ iBo/b£ r\j Lz as discussed in Section EH 
















14 


Rappazzo 




Figure 7. Runs C (m=l) and D (m=0-4): Total energy vs. time for line-tied simulations with different values of feo/^o for run C with 
single parallel mode m=0 {left panel) and run D with m=l {right panel). The 2D run A with bo/Bo = 10% is added for comparison. The 
insets show in logarithmic scale total energies normalized with their initial values. 


Initially also the 3D line-tied system starts to behave 
as in the 2D case, with the tension of perpendicular field 
lines creating an orthogonal velocity, that coupled with 
all others nonlinear terms are the only ones that can cas¬ 
cade energy and generate current sheets. But this dis¬ 
places the total line-tied (axially directed) field lines, that 
now cannot be freely convected around as in the periodic 
case because of the line-tying constraint at the bound¬ 
aries, and is then counteracted by the enhanced axial 
tension that resists bending. Magnetic fields with inten¬ 
sities smaller than the threshold (|23|) a small variation of 
b along 2: (corresponding to a variation scale larger than 
the loop length Lz) is enough to reach an equilibrium, 
but for intensities larger than of 6 > iB^jLz current 
sheets must form and energy dissipate in order to reach 
the physically accessible equilibria with b ^ iBo/Lz^ 
since for larger magnetic field intensities the equilibria 
are strongly asymmetric along ^ and therefore physically 
inaccessible. 

The analysis in this section has considered exclusively 
initial conditions invariant along the ^-direction, with 
only the parallel mode m=0. In the following sections we 
extend it to include field variations along 2: with higher 
parallel modes.. 

4.2. Runs C: single mode m=l 

The finite length of coronal loops renders the system 
akin to a resonant eavity. A forcing velocity with fre¬ 
quency n at the photospheric boundary, e.g., 

u{x,y,z = L,t) =Uph(x, y) cos{Lot) , (24) 

where uj = 27ru is the angular frequency, will inject 
Alfven waves at that frequency propagating in the ax¬ 
ial direction of the loop. In general these waves, that 
are continuo usly injected and reflected at the boundaries 
(see Section [3Tl Equation [IS]), will be out of phase and 
decorrelated along the loop so that their sum will re¬ 
main limited the whole time to values of the order of the 
forcing velocity at the boundary, with no growth in time 
for the amplitude of the resulting velocity and magnetic 
fields. But for the resonant frequeneies 

z^n = ^^A/2, with n G N (25) 


and ua = I/ta, the wav es are in phase and they sum 
coherently (|Ionsonl ri985[) . Thus the magnetic and ve¬ 
locity fields in the loop grow linearly in time similarly 
to the case with constant velocity (Equation [IS]). Eor 
instance, considering the boundary velocity (|M|) at the 
resonant frequency with n > I, the resu l ting fields 
grow approximately as (lEinandi fc Velli]ll999l: iRaopaz^ 
l2(M IChinderi VellilDfUlb 

b ^ Uph cos ( CJn — ) COs{uJnt) — , (26) 

\ VaJ ta 

u ^ Up/, sin [uJn—] sm{uJnt) —. (27) 

\ VaJ ta 

Indeed a constant velocity can be regarded as the zero 
frequeney resonanee n = 0 of the system, that differs 
from resonances with n > I because while the magnetic 
field grows linearly in time, for n = 0 the velocity field 
does not grow and its value remains of the same order of 
magnitude of the photosphe ric v elocity {u ^ Uph). 

As mentioned in Section EH for X-ray bright solar 
coronal loops photospheric motions have a low frequency, 
giving rise to a coronal magnetic field dominated by the 
parallel m = 0 mode, and their low frequency can then 
be approximated with zero. In general photospheric mo¬ 
tions characterized by a surface convective timescale Tsc 
will not have a single harmonic at the frequency Ijrsc^ 
rather the amplitude of its Eourier transform will peak at 
the frequency 1/tsc but include many other harmonics. 

Eor longer loops (with longer Alfven crossing times 
Ta = LzjvA comparable or longer than granulation 
timescales), and for loops on other active stars with mag¬ 
netized coronae and outer convective envelopes, photo¬ 
spheric motions can have frequencies closer to resonances 
higher than z^o = 0- Eurthermore also when photospheric 
motions have a dominant low frequency, higher frequency 
modes will be present, although with smaller ampli¬ 
tudes (iNigro et al.l l2QQ8f) that contribute considerably 
less heating (jMilano et al.l 1 199 71 ). In all these cases pho¬ 
tospheric motions will give rise to parallel modes higher 
than zero (m > I) for the coronal magnetic field. These 
will be the dominant modes when the photospheric fre- 
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Figure 8. Runs C (m=l) and D (m=0-4): Magnetic field lines of the orthogonal magnetic field component b and current density j in the 
mid-plane 2 ; = 5 at selected times for runs C (top row) and D {bottom row) with bo/Bo = 10%. The left panels show the initial condition 
at t = 0, the central panels show the fields at the time of maximum dissipation, while the right panels show the fields at a later time when 
they have relaxed (asymptotic regime). 


quency is resonant, and give a small contribution to the 
magnetic field when photospheric motions frequencies are 
close to zero. Additionally higher parallel modes can also 
be generated by nonlinear dynamics also when st arting 
with a zero parallel mode (|Bnchlin fc Vellil I2QQ7I) . and 
by disturbances stemmi ng from chromospheric dynamics 
(|De Pontien et al.ll2QQ7[) . 

It is therefore of interest to consider initial conditions 
with modes higher than m = 0, and in the numerical 
simulations analyzed in this section the initial magnetic 
field has only the parallel mode m = 1 (corresponding 
to the resonant frequency 1^2 = while large-scale 

perpendicular modes are set as in previous simu lations 
with wavenumbers between 3 and 4 (Section 12.11 Equa¬ 
tion 0 )- Thus unlike run B with m = 0, the out-of- 
equilibrium initial magnetic field now varies also along 
with both terms in the equilibrium equation ([8]) not 
vanishing {dzj 7^ 0, b • Vj 7^ 0). The magnetic field 
lines of the orthogonal component b and current den¬ 
sity j are shown in Figure [8] (top row) at time t = 0 
in the mid-plane z = 5, while isosurfaces of the mag¬ 
netic potential at t = 0 are shown in Figure [9] (central 
column). In both figures the case with Bq/Bq = 10% is 
considered. This is one of a series of simulations collec¬ 
tively labeled as runs C, with same parameters for the 
initial condition except the magnetic field intensity bo 
(the multiplicative factor in Equation [3]) that spans the 
range 0.1% < bo/Bo < 10%. 

The time evolution of total energy for runs C with 
different values of bo is shown in Figured (left panel). 
The run with bo/Bo = 10% has a similar behavior to 
run B with same magnetic field intensity (cf. Figure [2j). 
Its energy decays approximately as F oc with current 
sheets forming in physical space (Figure [HI top row) and 


dissipating ^ 92% of the initial energy, a slightly higher 
value respect to the corresponding run B. Subsequently 
the system relaxes to an asymptotic state with b/Bo ^ 
2.87%. 

The analysis of the equilibria set forth in Section 13.11 
shows that the only dynamieally aeeessible equilibria are 
those with variation length-seale greater than approxi¬ 
mately the loop length zi > These have structures 
very elongated in the axial direction, and therefore dom¬ 
inated by the parallel mode m = 0. When initial condi¬ 
tions do not include the m = 0 mode, their higher modes 
will necessarily have to transfer part of their energy to 
the mode m = 0 via nonlinear dynamics in order to relax 
to equilibrium. Furthermore most of the energy of the 
modes with m > 1 that is not converted into the parallel 
zero mode must be dissipated for the relaxed state to be 
close to a reduced MHD equilibrium with Z£ > Lz, with 
a predominant parallel zero mode. 

In fact the isosurfaces of the magnetic potential tp in 
Figure|9]show that both runs B and C with bo/Bo = 10% 
(respectively in the left and central columns) relax to a 
lower energy state with structures very elongated along 
z (the computational box has been rescaled for an im¬ 
proved visualization, but the axial length is ten times 
longer that the perpendicular cross section length), even 
though their initial conditions are radically different, con¬ 
sisting exclusively of mode m = 0 for run B and m = 1 
for run C. 

Consequently similar dynamics will occur for all runs C 
independently from the value of bo/Bo, as shown in Fig¬ 
ure [71 because all of them do not have a mode m = 0 
in their initial magnetic field. Thus they all decay while 
part of the energy initially in the m = 1 mode is either 
transferred to the m = 0 mode (and partially also to 
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higher order modes) or dissipated (including dissipation 
of the higher order modes generated during the nonlinear 
dynamics). On the contrary when the initial condition 
is made only of mode m = 0, the system is very close to 
an equilibrium for bo/Bo < 3% (runs B, Figure [2]) since 
this condition implies Z£ > Lz and the magnetic field is 
quasi-invariant along z, so that no substantial dynamics 
occur when Bq/Bq < 3%. 

Figured shows that for runs C energy starts decaying 
at later times for smaller values of bo/Bo. The longer 
timescales for dissipation to occur and energy to start 
decaying is consistent with the decrease of the strength 
of nonlinear interactions for lower values of bo/Bo. For 
instance the eddy turnover time increases as ti ^ 21/bo, 
since no velocity is initially present and this soon be¬ 
comes of the order of bo/2 (because for resonant fre¬ 
quencies velocity is in equipartition with the magnetic 
field). Subsequently velocity strongly decreases once a 
zero mode is created and the system relaxes. 

Clearly the fraction of magnetic energy dissipated dur¬ 
ing the decay depends on several factors. The most im¬ 
portant of these is how much energy in transferred to the 
parallel zero mode, since all higher modes will be largely 
dissipated in order to reach an equilibrium with zi ^ Lz. 
A detailed analysis of the energy fluxes between these 
modes goes beyond the scope of the present paper, and 
might be carried out in future work. Nevertheless it is 
clear from FigureElthat runs C with 4% < bo/Bo < 10% 
generate a zero mode quickly, while for bo/Bo < 1% a 
large fraction of the initial energy is dissipated with only 
a small fraction transferred to the zero mode. In all cases 
in the asymptotic stage, when the mode m = 0 is the 
strongest mode, the relaxed magnetic field intensity b is 
below the stability threshold (|2^ with b/Bo ^ 3%. 

In spite of all the aforementioned differences, the longer 
decay timescales for runs C with lower values of bo/Bo 
render similar the behavior of the system forced by 
photospheric motions for both cases with a velocity 
that is constant in time (zero frequency) and a veloc¬ 
ity with higher resonant frequencies. In fact considering 
a straightened loop with initially only the guide field Bo, 
if a constant velocity (pq = 0) is applied at the photo- 
spheric boundaries, the magnetic field will grow linearly 
in time initially (Equation [E]), because until the or¬ 
thogonal component of the magnetic field does not reach 
the critical value b ^ iBo/Lz (Equation [23]) the sys¬ 
tem is very close to equilibrium and nonlinear terms can 
be neglected. The linear growth of the magnetic field is 
derived indeed from the linearized reduced MHD equa¬ 
tions ([T^ . 

In similar fashion if the photospheric velocity frequency 
is a higher resonance Vn = '^^a/2, with n > I, again non¬ 
linear terms can be neglected initially because for low val¬ 
ues of the magnetic field intensity b the decay timescales 
are much longer than the linear growth of the mag¬ 
netic field, with the amplitude doubling every Alfvenic 
crossing time ta. Equations are obtained also 

for the resonant frequencies from the linearized reduced 
MHD equations, analogously to the constant forcing case 
(Equations [I5]-[T6]). A statistical steady state will fi¬ 
nally be obtained when the energy flux injected in the 
system at the boundary by photospheric motions is bal¬ 
anced by a similar energy fiux from the large toward the 
small scales (to form current sheets), in similar fashion 


to the constant velocity case (|RaDDazzo et al.|[2QQ7l) . 

4.3. Runs D: modes m=0-4 

In this section we analyze numerical simulations 
(runs D) in which the initial magnetic field b includes all 
parallel modes me [0,4], while large-scale perpendicular 
modes are set as in previous si mula tions with wavenum¬ 
bers between 3 and 4 fSection 12. li Equation [3]). The 
parallel zero mode is the strongest, with higher parallel 
modes having less energy. Eor all runs the fraction of 
magnetic energy in the parallel mode m, indicated with 
(om, is set in the initial magnetic field (Equation |3]) so 
that ^rn/<^o = (^ + I)~^‘^5 corresponding to a progres¬ 
sively smaller energy for higher modes. Explicitly 
= 16.5%, 5.7%, 2.7%, 1.5% for m e [1,4]. This spe¬ 
cific choice of values is arbitrary, but the presence of 
multiple parallel modes of decreasing weight is chosen 
to represent the coronal field induced by photospheric 
motions in which multiple frequencies are present. Al¬ 
though to date there are no measurements of the spec¬ 
trum of photospheric velocities, the presence of higher 
frequency modes with decreasingly smaller amplitudes is 
expected. Since smaller granules are expected to have 
faster convective timescales this is partially confirmed 
by the recent detection of mini-granular structures with 
their size distributed as a power law with an approxi¬ 
mate Kolmogorov in dex ^ — 5/3 and dominan t on scales 
smaller than 600 km (| Abramenko et al]l20I2[) . 

As shown in Eigure [7] (right panel) the dynamics are 
analogous to those of run B with only the parallel mode 
m = 0 (cf. Eigure [2|). The dissipated energy decreases 
for lower values of bo/Bo, but unlike runs B a small but 
discernible energy dissipation occurs also for very small 
ratios of bo/Bo. Each run D dissipates a larger fraction 
of magnetic energy respect to the corresponding run B 
with same bo /Bo before reaching the asymptotic regime 
(cf. insets in Eigures[7]and[2|), because initially a fraction 
of the energy in runs D is in parallel modes higher than 
zero, and they will be dissipated during the decay so that 
the system can relax to equilibrium (dominated by the 
mode m = 0 with variation length-scale Z£ > Lz). 

Dynamics in physical space are also similar to those of 
run B, as shown in Eigure [8] (bottom row) for bo/Bo = 
10%, with current sheets forming and dissipating en¬ 
ergy thus leading to a relaxed field with larger magnetic 
islands through an inverse cascade. The evolution of 
the three-dimensional structure of the magnetic poten¬ 
tial shown in Eigure [9] (right column), is similar to 
those of runs B and C. The relaxed magnetic potential 
at t = 240 Ta has a structure elongated along z, with a 
strong m = 0 parallel mode similarly to runs B and C, 
even if its structure at time t = 0 is considerably more 
complex due to the presence of multiple parallel scales 
(m e [0,4]). 

5. CONCLUSIONS AND DISCUSSION 

Equilibria and dynamics of the magnetically confined 
regions of solar and stellar atmospheres have been inves¬ 
tigated with a reduced MHD cartesian model to advance 
our understanding of the mechanism that powers the X- 
ray activity of the Sun, late-type main sequence stars, 
and more in general of stars with a magnetized corona 
and an outer convective envelope. 
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Figure 9. Isosurfaces of the magnetic potential 'll; (in yellow and transparent red and blue colors) for the three-dimensional simulations 
with line-tied boundary conditions runs B, C and D, with respectively single parallel modes m = 0, m = 1 and all modes m G [0,4], and 
initial conditions with 6o/^o = 10%. The elongated structure of 'ijj along 2 : shows that in the asymptotic regime the fields relax into an 
equilibrium with a strong parallel m=0 mode not only for run B, that initially has only the mode m=0, but for all the initial conditions 
considered, including runs C that initially has only the parallel mode m=l and run D that is started with many parallel modes. Snapshots 
times, left to right, are respectively t = 516 ta, 893 ta and 240 ta- The magnetic potential is defined in Section [2] The computational 
box has been rescaled for an improved visualization, but the axial length is ten times longer that the perpendicular cross section length. 


Since equilibria play a pivotal role in understanding 
the dynamics of this system, their structure has been 
analyzed in detail in Section [3l The mapping between 
the solutions of the 2D Euler equation U2i:)(x, t) and re¬ 
duced MHD equilibria beg(x^,z) (Equation [TO]: t ^ 
U 2 d beg/^o) allows to formulate a heuristic quantita¬ 
tive analysis of the structure of the reduced MHD equi¬ 
libria. The inverse easeade developed by the solutions of 
the Euler equation in time, corresponds to an asymmet- 
rie strueture along z of the reduced MHD equilibria, as 
pictorially summarized in EigurelH 
In similar fashion to 2D velocity vortices of scale i that 
double their size in about one eddy turnover time t£ ^ 
^|u£ (Equation [l3]), a reduced MHD equilibrium with 
orthogonal magnetic field of intensity b^d at the bottom 
boundary z = 0, and made of magnetic islands of scale £, 
will have progressively larger magnetic islands at larger 
z, doubling their transverse scale over the axial spatial 
distance 

(28) 

Obd 

This represents the parallel variation length-seale of 
the equilibrium solution, and measures quantitatively its 
asymmetry. An equilibrium is strongly asymmetric along 
z if the variation scale is smaller than the loop length 


Z£ < Lz, but if the variation scale is greater than approx¬ 
imately the loop length Z£ > Lz then the equilibrium is 
quasi-invariant along 2; (at the loop scale). 

On the other hand in reduced MHD any spatial vari¬ 
ation of the physical fields along 2; is rapidly propagated 
away at the Alfven speed Bq (the fastest speed in the 
system) by the linear terms (x BqOz in Equations ([I])-(j2j). 
Therefore the physieal solutions eannot develop strongly 
asymmetrie struetures along 2(, as confirmed by bound¬ 
ary forced and decaying numerical simulations with 
line-t ying boundary conditions and a strong guide field 
(e.g.. Galsgaard &: Nordlundl 119961: iDmitrnk fc Gomezl 
I999I: ^DPazzo et al.l I2QQ8I: IWilmot-Smith et al.l I2Q11I: 
Dahlbnrg et al.ll2Q12[) . This implies that physical solu¬ 
tions are close to an equilibrium or not depending on the 
intensity of the orthogonal component of the magnetic 
field 6, that regulates the value of the parallel equilibrium 
variation scale Z£ ([28|) for a given guide field of intensity 
Bq. Gonsequently the redueed MHD equilibria that ean 
be aeeessed dynamieally are those quasi-invariant along z 
with an orthogonal magnetie field eomponent b for whieh 
the parallel length-seale is greater than approximately the 
loop length Z£ ^ Lz. 

The simulations reported by iRappazzo fc Parker! 
(HoH), analyzed in detail in Section 14.11 (runs A-B), 
are in agreement with this picture. They considered ini- 
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tial magnetic fields bo invariant along 2; (i.e., only the 
parallel mode m = 0 is present initially and djo = 0), 
but with different intensities and non-vanishing Lorentz 
force for the perpendicular magnetic field, i.e., bo • Vjo 7^ 
0, and identified the magnetic field intensity threshold 
b iBojLz, corresponding to the critical variation scale 
Zi ^ Lz. Initial conditions with bo > iBojLz (z£ < Lz) 
have dynamics increasingly similar to 2D MHD turbu¬ 
lence decay for larger values of bo^ with the orthogonal 
magnetic field forming current sheets and dissipating en¬ 
ergy that decays as a power-law in time with E oc t~^ 
{a 1 for bo/Bo = 10%, Figure [2]). As shown in 
Equation ([18)) in the 2D case the evolution of fields 
with different initial intensities is self-similar in time 
with b'(t) = b(t • 60/60)60/60, implying that they all 
decay with same power-law index and relax to asymp¬ 
totic fields with intensities proportional to their initial 
value b'^/b'o = boo/bo. But in stark contrast with the 2D 
case the line-tied 3D simulations decay with progressively 
shallower power-laws (Figure [2]) for weaker initial mag¬ 
netic fields with smaller ratios bo/Bo, relaxing to asymp¬ 
totic equilibria for which the ratio boo/bo is not indepen¬ 
dent from bo. Instead boo ^ ^Bo/Lz, corresponding to 
a variation scale approximately equal to the loop length 
Z£ ^ Lz (as explained in Section H7T] the orthogonal scale 
i increases for stronger magnetic field because an inverse 
cascade of magnetic energy occurs). Furthermore ini¬ 
tial magnetic fields with intensity below the threshold 
bo ^ £Bo/Lz {z£ > Lz) show little dynamics with no 
significant decay nor current sheets formation and dissi¬ 
pation. 

Therefore these simulations confirm numerically that 
the dynamically accessible equilibria are those quasi¬ 
invariant along 2; (i.e., with a dominant m = 0 mode) 
with magnetic field intensity smaller than the thresh¬ 
old 6 < iBo/Lz, and corresponding parallel variation 
length-scale larger than approximately the loop length 
Z£ > Lz. The nature of this equilibria is radically dif¬ 
ferent from the classic reduced MHD equilibria consid¬ 
ered in plasma and solar physics in the framework of lin¬ 
ear instabilities (kink, tearing, etc.), that typically are 
strictly invariant along 2; {dz = 0) and in the reduced 
MHD framework have a vanishing orthogonal Lorenz 
force component with b • Vj = 0. This condition is 
satisfied by very symmetric fields, e.g., a sheare d field, 
or circular field lines (examples can be found in iParkerl 
119831: iLongcope fc Stransslll993L and refe renc es therein). 
But our initial magnetic fields (Section 12.ip have non¬ 
vanishing orthogonal Lorentz forces (b • Vj % 0), a prop¬ 
erty that stems from the complexity and disorder of pho- 
tospheric motions. As shown in Figure [5] and [6] both 
terms in the equilibrium Equation ([8]) do not vanish as 
the system relaxes to equilibrium, with their rms az and 
(Equation [20]) getting asymptotically equal, while 
the rms of their sum aeq (Equation [2T]) vanishes asymp¬ 
totically as a power-law. Consequently the system does 
not relax to a classic linearly unstable equilibrium with 
dzj = 0 and b • Vj = 0, as confirmed by a visual inspec¬ 
tion of the orthogonal magnetic field b in Figure |4] (right 
column). 

The simulations (runs B) have very different dynamics 
whether their initial parallel variation scales are larger 
or smaller than the critical length-scale Z£ ^ Lz. For 
Z£ < Lz, with 6 > iBo/Lz, the initial magnetic field 


is very far from the corresponding equilibrium, as this 
is too asymmetric along ^ and is therefore dynamically 
inaccessible. The only way for the out-of-equilibrium 
field to reach an equilibrium is therefore to decay to a 
lower energy configuration with smaller 6 until the criti¬ 
cal length scale Z£ ^ Lz is reached, and this necessarily 
implies the formation of current sheets and dissipation 
through nonlinear dynamics (a magnetically dominated 
nonlinear MHD turbule nt cascade analyzed in depth in 
iRappazzo fc Vellil l2Qll[) . On the contrary initial mag¬ 
netic fields with Z£ > Lz, for which 6 < iBo/Lz, are 
very close to the corresponding equilibrium because they 
are both elongated along z. Thus the field simply read¬ 
justs to the close equilibrium with no significant nonlin¬ 
ear dynamics, current sheet formation nor dissipation, as 
shown in Figures EHD and particularly in Figures [5] and 
|6l Strictly speaking also in this case a very small energy 
dissipation occurs, but it is negligible, does not involve 
the formation of significantly stronger currents, and addi¬ 
tionally nonlinearities are diminished in close proximity 
to equilibrium. 

The quasi-static evolution of the magnetic field is then 
restricted only to field intensities smaller than approxi¬ 
mately 6 < iBo/Lz, while stronger fields are necessarily 
out-of-equilibrium and develop turbulent dynamics with 
subsequent current sheet formation and energy dissipa¬ 
tion. 

Consequently two distinct stages can be identified in 
the dynamics of an initially uniform and strong axial 
magnetic field Boez shuffled at its footpoints by a con¬ 
stant or low frequency photospheric velocity Uph (see Sec- 
tion l3.1l for a discussion on forcing frequencies). To mimic 
the solenoidal component of the photospheric horizon¬ 
tal velocity (the irrotational component cannot twist the 
field lines), the incompressible velocity at the boundary 
Unh i s made up of distorted vortices (see IRappazzo et M] 
I2QQ8L for a specific example) with Uph • ^ojph % 0, given 
the general complexity and disorder of photospheric mo¬ 
tions. At first photospheric motions generate an orthogo¬ 
nal coronal magnetic field component that grows linearly 
in time and is a mapping of the photospheric velocity, i.e., 
b = Uph t/xA (Equation [H]). Until its intensity remains 
below the threshold 6 ^ iBo/Lz the field is essentially 
in equilibrium and nonlinearities do not develop, leading 
to its linear growth in time. In fact neglecting nonlinear 
terms in the reduced MHD equations, the linear growth 
follows from the remaining linear terms (Equation [T5] ) 
and boundary conditions. When the magnetic field in¬ 
tensity crosses the threshold the variation length-scale of 
the corresponding equilibrium becomes smaller than the 
loop length Z£ < Lz. The structure of the equilibrium 
becomes then too asymmetric along z, while the dynam¬ 
ically induced magnetic field (Equation [16]) is quasi¬ 
invariant along z. The magnetic field is then too distant 
from its corresponding equilibrium that cannot be ac¬ 
cessed dynamically unless the field intensity decreases. 
The Lorentz force components that were in equilibrium 
during the quasi-static stage now cannot reach a balance 
with each other, the magnetic field is therefore in non- 
e quilibrin r n and nonlinear dynamics develop. 

(Parkerl ()1988f ) had conjectured a two-stage process 
to account for the inferred Poynting finx in active re¬ 
gions, estimated by IWithbroe fc Novesl ()1977[ ) at about 
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Sz ^ 10^ erg cm“^ s“^. In fact if current sheet for¬ 
mation and energy dissipation would be effective at an 
earlier stage, with too weak magnetic fields, then the 
Poynting flux would be too small to sustain the X-ray 
activity of active regions. Reverting now to conven¬ 
tional units, the time and space averaged Poynting flux is 
given by {Sz) ^ S/P ^ P^a,\\ (see Section [221 

density p and Alfven velocities are introduced from di¬ 
mensional calculations), where va^\\ = Bq/ pAnp and 
va,i. = h/PAnp are the Alfven velocities associated re¬ 
spectively to the guide field Bq and the orthogonal mag¬ 
netic field component b. Introducing the threshold mag¬ 
netic field intensity b ^ IBq /Lz (or the associated Alfven 
velocity) we obtain for the Poynting flux: 


{S^) ~ pv\ ^^Uph- 


Bq Up}i£ 

47rI/2 


(29) 


This coincides with the strong guide field reg i me o f 
the scaling relation obtained by iRaPDazzo et al.l (|2QQ8f ) 
(Equation [68] with a ^ I) for boundary forced simula¬ 
tions, that yields for typical solar active region loops an 
energy flux ^ 1.6 x 10^ erg c m~^ s~^, in the lower rang e 
of the constraint inferred by IWithbroe fc Nove'sl (119771 ) . 
But recent fully compressible MHD simulations with sim¬ 
ilar setup, that include the integration of an energy equa¬ 
tion with thermal conduction and energy losses provided 
by optically thin radiation, and in addition have den¬ 
sity stratification (with strong gradients from the chro¬ 
mosphere to the corona), exhibit Poynting fluxes of the 
order of ^ 10^ erg cm“^ s“^, and most importantly an X- 
ray emission that matches the physical properties of the 
observed radiation (jPahlburg et al.ll2QI5. submitted^ . 

Additionally numerical simulations with initial mag¬ 
netic fields not invariant along z have been carried out 
(runs C and D). Higher parallel modes can be present for 
both the Sun and other active stars of interest. They can 
be generated by the nonlinear dynamics even when pho- 
tospheric motions have a low frequency, or they can be 
directly excited by photospheric motions in long loops, 
that on some stars have observation ally inferred lengths 
of the order of several stellar radii (iFavata et al.l I2QQ5I: 
iGetman et al.ll2QQ^ [Peterson et al.ll^IQl) . 

Runs C include only the parallel mode m = I, while 
for runs D all modes me [0,4] are present. As in previ¬ 
ously discussed runs B the initial magnetic field bo is not 
in equilibrium, now with both terms in the equilibrium 
Equation ([8j) non-vanishing (bo • Vjo 7^ 0 and dzjo 7^ 0). 
Remarkably these initial conditions decay to equilibria 
with structures similar to those of runs B (whose initial 
magnetic fields included only the mode m = 0), as shown 
in Figure [9l Independently from the structure of the ini¬ 
tial magnetic fields and from the specific modes that it 
includes, the final equilibrium is always quasi-invariant 
along z with parallel variation length-scale larger than 
approximately the loop length zi > Lz. Its structure 
is elongated along the axial direction, a strong m = 0 
parallel mode is present, and the magnetic field inten¬ 
sity is smaller than the threshold value b < IBq/Lz. 
This further confirms the analysis of the equilibria per¬ 
formed in Section ixn i.e., that the reduced MHD equi¬ 
libria dynamically accessible in a line-tied configuration 
are elongated along the axial direction with a domi¬ 
nant m = 0 parallel mode. It also implies that non¬ 


linear dynamics can transfer energy from higher parallel 
modes to the mode m = 0 even when this is not initially 
present (runs C), a p rocess that can be of interest also in 
perio dic turbulence (|Alexakis[ [2QII[: [Schekochihin et ahl 
[2QI2E and has been conjectured to play a role in the 
dynamics that lead t o the acceleration of the solar wind 
(|Dmitrnk et al.[[2QQI[) . 

In contrast to runs B, with initial conditions invari¬ 
ant along 2), that decay only for magnetic field intensi¬ 
ties above the threshold b ^ IBq/Lz, in runs C a decay 
is always observed independently from the intensity of 
the initial magnetic field (Figure [7j). Since the accessi¬ 
ble equilibria are quasi-invariant along 2: and the initial 
magnetic field of runs C does not contain the m = 0 
mode but only the m = 1 mode, then it always decays. 
The asymptotic stage is reached when the mode m = 0 
has been generated and excess energy in higher modes 
dissipated. The intensity of the relaxed magnetic field 
depends on the ability of nonlinear dynamics to trans¬ 
fer energy among higher parallel modes and from these 
to the m = 0 mode, but the longer decay timescales 
for lower values of bo/Bo are consistent with a decrease 
of the strength of nonlinear interactions (e.g., the eddy 
turnover time decreases ad ti ^ ^i/bo^ see Section [T2]) . 
The longer nonlinear timescales render the effect of a 
high-frequency resonant photospheric forcing similar in 
many aspects t o th at of a constant photospheric veloc¬ 
ity (see Section [T2] for a more complete discussion). In 
fact if a photospheric velocity with a higher resonant fre¬ 
quency Un = nvAl‘^ with n > I is applied at the bound¬ 
ary, nonlinear terms can be neglected initially because 
for low values of the magnetic field intensity b the de¬ 
cay timescales are much longer than the linear growth 
of the magnetic field. A statistical steady state will fi¬ 
nally be obtained when the energy flux injected in the 
system at the boundary by photospheric motions is bal¬ 
anced by a similar energy flux from the large toward the 
small scales (to form current sheets), in similar fas hion 
to the constant velocity case (|RaDDazzo et al.[[2QQ7l) . 

When initial conditions include higher parallel modes 
m G [0,4] and the m = 0 mode has the largest amplitude 
(runs D, Section 14. 3p the dynamics are very similar to 
runs B with only the m = 0 mode in the initial conditions 
(cf. Figures[2]and[7|). The excess energy in higher parallel 
modes is either dissipated or transferred to the m = 0 
mode and the relaxed fields have structures similar to 
runs B and C (Figure [9]). 

The parallel variation scale zi ^ iBo/b introduced here 
(Equation [28]) can be interpreted as a critical length or 
twist. In fact given the magnetic field intensity 6, sig¬ 
nificant nonlinear dynamics will develop only if the loop 
length is longer than the variation scale Lz On the 


other hand, fixed the loop length Lz nonlinear dynamics 
will develop only if the variation scale is smaller than the 
loop length Z£ < Lz^ or equivalently if the field intensity 
is larger than the threshold b > IBq/Lz^ that corresponds 
to an average twist larger than (4>) ^ Lzb/{iBo) > n/3 
(this is only an estimate, since the orthogonal field does 
not have cylindrical symmetry b 7^ b{r), the twist should 
be computed numerically for sample field lines). 

The concept of a critical length or twist has been 
developed in the study of several linear instabilities 
in coronal loops with line-tied boundar y cond i tions, 
including kink and tearing instabilities ([Raadd [I972[: 
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Hood & Priest! [1979L [1981!: [Einaudi & van Hoven[[198lL 

198^ 

UVelli & Hood! 1989!: [Velli et al.!! 1990!: [Foote ^ 

^ Craid 


[ Lionello et ah! [1998!: [Huane: & Zweibel! 

to 

0 

0 


Huang et al.ll2QlQ[ r It is found that the system is lin¬ 


early unstable for a fixed magnetic field intensity only if 
the loop length is larger than a critical value (for kink 
and other instabilities a critical twist can be used equiv¬ 
alently) . 

It is important to remark that although it is useful 
to regard the parallel variation scale Z£ (Equation [28] ) 
as a critical length, and that an average critical twist 
can be defined, the dynamics that they help describe 
are not linear instabilities. In the configurations of in¬ 
terest to this paper the critical length or twist distin¬ 
guish two different dynamic regimes in which respectively 
for small field intensities b nonlinear dynamics are sup¬ 
pressed and the evolution can be regarded as quasi-static, 
while for stronger intensities nonlinear dynamics develop. 
Of course the boundary between these two regimes is not 
sharp, but gradually as b is inereased the Lorentz 

foree b • Vb eannot be balaneed by the axial field line ten¬ 
sion B^dzh, leading the system out of equilibrium. In 
particular all the 3D line-tied magnetic fields that we 
have considered lack the symmetries of linearly unstable 
configurations, and for all of them the orthogonal com¬ 
ponent of the magnetic field is not symmetric and its 2D 
Lorentz force component (for which b • Vj 7 ^ 0) does not 
vanish at all times, from the initial condition to the re¬ 
laxed asymptotic stage, as can be seen in Figures |4][6l [HI 
andlH 

The reason for which the equilibria with Z£ > Lz are 


not linearly unstable is that they are close to each other. 
Therefore adding a perturbation to the magnetic field 
simply changes slightly the corresponding equilibrium to 
which the field readjusts. For instance the initial condi¬ 
tion of run B with bo/Bo = 2% is not exactly an equilib¬ 
rium (since dzjo = 0 and bo • Vjo 7 ^ 0 ), therefore it can 
be regarded as an equilibrium to which has been added 
a small perturbation. But as shown particularly well 
in Figure [ 8 ] (bottom panel) no instability of sort is de¬ 
tected, rather the field undergoes a slight readjustment. 
Clearly for a non-vanishing orthogonal magnetic field b 
in equilibrium, with a progressively smaller 2D Lorentz 
force component for which b • Vj ^ 0 (with 6 7 ^ 0 ), 
the field approaches a symmetric configuration that can 
be linearly unstable, since from the equilibrium Equa¬ 
tion dH]) also dzj ^ 0 in this limit. But as previously dis¬ 
cussed the configurations of interest here are those with 
non-vanishing orthogonal 2D Lorentz force component 

b • Vj 7 ^ 0 . 

The relaxation of coronal magnetic fields has of¬ 
ten been studied subsequently to a linear instabil¬ 
ity, mostly kink modes. Particularly for the cases 
that have a strong axial magnetic field the structure 
of the lo wer energy relax e d field appears elongated 
along 2 ; (iMikic et ah! 1 19901: iLongcope fc Strauss! 1994 


Baty fc Hevyaertsll19961: IVelli et cd 


Batv 2QQQI: iGerrard et al.l I2QQ2I: [Browning et ah. 
Hood et ni [2QQ9[) . Early b oundary forced sim- 


Lionello et al 


1998 : 

2008; 


ulations have been performed by [Ng fc Bhattachar'i^ 

S , with similar setup as those of [Rappazzo et al. 

. But due to their low resolution they misinterpret 
as an instability the dynamics that develop as the thresh¬ 
old b ^ IBojLz is crossed, when the system gradually 


transitions from quasi-static evolution to turbulent dy¬ 
namics. From the simulations (runs B) and the equilibria 
analysis it is clear that the forces that are approxima- 
tively in balance below the threshold become gradually 
unbalanced for larger magnetic field intensities, leading 
to the development of nonlinear dynamies with no in¬ 
termediate instability as would oeeur for instanee with a 
kink mode, where the nonlinear stage would follow the 
linear instability. The non-vanishing 2D Lorentz force 
term b • Vb (for which b • Vj 7 ^ 0), that in the 2 D 
case (run A) sets the system out-of-equilibrium and de¬ 
velops turbulent nonlinear dynamics, can be balanced 
in the 3D line-tied runs B by the axial field line ten¬ 
sion term Bodzh for field intensities below the threshold 
b < iBojLz. For larger magnetic field intensities the 
2D Lorentz force term is stronger than its axial com¬ 
ponent, hence the dynamics develop progressively more 
akin to the 2 D case. Ultimately a force balance cannot 
be reached because the corresponding equilibrium is too 
asymmetric along 2 : and therefore dynamieally inaeees- 
sible, so that for larger intensities b a larger fraction of 
energy must be neeessarily dissipated for the system to 
be able to relax and access a new equilibrium. 

This two-stage process then provides a fully self- 
consistent alternative to coronal heating models based 
on instabilities. For instance, since resistive instabili- 
ties are slow for macroscopic ally thick magnetic shears, 
[Pahlburg et al.[ ([20051 [20091 ) obtain a shear intensity 
threshold for dissipation supposing that nanoflares would 
occur when photospheric motions shear the magnetic 
field beyond a certain angle, when a secondary ideal 
instability (triggered by the slow primary resistive in¬ 
stability) can develop thus accelerating the dynamics. 
On the other hand, as discussed in this paper, as 
photospheric motions disorderly twist the field lines, 
once the magnetic field intensity is higher than the 
threshold b > IBojLz magnetically dominated turbu¬ 
lent dynamics develop, forming current sheets that thin 
down to the dissipative sca les on fast Alfven time-scales 
([Rappazzo fc Par]^[2Q13[ ). while tr iggering the “Meal” 
tearing instability ( see Introduction ; [Pucci fc Vellil[2Ql3: 
[Landi et~an [2015!: [Tenerani et al.[ [2Q15[ ). and leading 
to dynamics similar to so-called plasmoid instability 

([Bulanov et al.[[T978[:[BiskamD[[l98^lLoureiro et al.[[2QQ7i: 

[Lapental[2QQ8[: iBhattachariee et al']i2()Q9[). 

Finally, these simulations of decaying magnetic fields 
show that, beyond the intensity threshold [Equa¬ 
tion ([28]) ]. current sheets form on fast ideal timescales 
beeause of the nonlinear dynamies that develop. This is 
in stark contrast with the frequent hypothesis of quasi¬ 
static evolution of the coronal magnetic field subject to 
footpoint shuffling, that should continuously relax to a 
nearby equilibrium wit hout forming current sheets (e.g., 
[van Ballegooiien[[l985[ ). In the quasi-static scenario the 
corona could be heated by the uniformly distributed 
small-scale current sheets created by the shredding of 
the coronal magnetic field after many s uccessive random 
walk steps of its field lines footpoints ([van Ballegooiien! 
[1986[ ). But this mechanism would lead to current sheet 
formation on timescale longer than photospheric con¬ 
vection (several random walk steps would be required). 
While the relaxation simulations pre sented in this pa¬ 
per and in [Rappazzo fc Parker! ([2013! ) show that current 
sheets form on ideal Alfven timescales (much faster than 
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convective timescales), with the footpoints fixed at the 
photospheric plates where no motions are in place (and 
therefore no footpoint random walk occurs). 
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